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ABSTRACT 


Power electronics have advanced to the point that they can be considered for use 
in large, high power dc electrical distribution systems. The proposed Navy DC Zonal 
Electrical Distribution System (DC ZEDS) being designed for the new DD-21 utilizes 
a rectified ac generator output which is filtered and stepped to usable voltages by 
local dc-dc converters. One characteristic of the high-bandwidth converters is a 
negative input impedance, which when coupled with an LC input filter, can lead to 
system instabilities. This thesis examines various stability criteria for determining 
parameters of the dc link filter. Comparisons between a simplified system model, a 
model using subsystem impedances and an Advanced Continuous Simulation 
Language (ACSL) model of a reduced-order system are made. Simulations were 
conducted to verify the validity of the stability criteria. The ACSL model provides 
an extremely useful tool in evaluating the response of various system parameters to 
changes in design values. The design criteria examined in this thesis can be 
ultimately used to provide design specifications to future vendors. 
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I. 


PROJECT BACKGROUND 


Many recent issues have evolved which warrant the need and provide the 
motivation for a new way of thinking about electric power distribution on-board Navy 
ships. In recent years, the Navy has pushed for more use of Commercial-Off-The-Shelf 
(COTS) parts in their new ship designs and existing ship upgrades. Hand-in-hand with 
this trend is the updating and designing of systems using the latest and most advanced 
technology available on the world market. Budgets are currently limited, which is 
forcing the need for cost effective designs. As ships are getting more complex, the 
possibility of automated systems is growing. The move towards automation is also 
driven by the economics of streamlining the manning of future Navy ships. Naval 
architecture is experimenting with unique designs which require innovative ship system 
designs to allow for more flexibility in their arrangements. Currently, the Navy 
implements a radial distribution configuration for electric power on-board its vessels. A 
suggested alternate option is that of zonal distribution. 

The Direct Current Zonal Electrical Distribution System (DC ZEDS) is a 
component of the Integrated Power System (IPS) currently being proposed for use on the 
new DD-21 design. The basic concept for the IPS is to use a common set of prime 
movers to generate power for both propulsion and distribution. It is fundamentally 
different from the radial design in many respects. These differences provide numerous 
advantages as outlined below. 

A. RADIAL DISTRIBUTION 

The radial electric power distribution concept involves distributing 60 Hz 
alternating current (ac) power produced by a limited number of generators to loads 
throughout the ship in a pyramid type configuration. Branching occurs from 
switchboards which route power to power panels located in the general area of the loads 
that they supply. Loads which are considered non-vital receive power from a single 
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source. Loads which are considered vital for ship operation are made more reliable by 
powering them from two different sources through merging branches of the pyramid. 
Power panels typically have multiple switchboard inputs. As a result, isolating power to 
a local area in a ship is unnecessarily complex. Standardization of the electric power 
distribution system among all ship classes is difficult at best due to the different number 
of power sources and vital loads. 

B. ZONAL DISTRIBUTION 

The zonal electric power distribution concept implements two power buses that 
run lengthwise through the ship in a configuration that, like the shipboard firemain 
system, maximizes reliability and minimizes the possibility that both buses are affected 
by a single battle casualty. All operating generators are connected to both buses via 
switchboards. Individual zones are designated typically to be areas which directly relate 
to flooding control zones. Localized power isolation would, therefore, be very clear-cut 
and easily learned during a sailor’s shipboard indoctrination period. Operation of a zonal 
distribution system would also be similar between ship classes. 

C. DC ZONAL DISTRIBUTION 

DC ZEDS expands upon the zonal distribution concept in the following manner. 
High voltage (approximately 1100 V) dc is distributed along the main power busses. The 
Ship’s Service Converter Modules (SSCM) are dc-dc converters which buffer the main 
bus and interzone loads as well as regulate output voltages. Where ac power is required, 
Ship’s Service Inverter Modules (SSIM) are used to convert the dc to ac and regulate the 
ac voltage and frequency within the zone. Vital loads are connected to both the port and 
starboard SSCM via a diode auctioneering process. Additional dc-dc converters 
synthesize the regulated lower dc voltages, for example 155 V and 270 V, required by 
combat systems. 
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The advantages of a DC ZED system are numerous. The most significant 
advantage is increased survivability during casualty situations. SSCMs can be 
programmed to act as “smart” circuit breakers. This would result in maximizing the 
ship’s capabilities to continue operating available equipment. A study has shown that a 
second improvement would be weight and space savings from the use of solid-state 
components and a reduced amount of cabling throughout the ship [1]. The use of a dc 
system simplifies operations such as paralleling loads. Only voltages need to be matched 
on a dc system, whereas voltages, frequency, and phase need to be matched on an ac 
system. A fourth benefit is that variable speed ac motor control is easily realizable with 
the Ship’s Service Inverter Modules (SSIM), which are programmable solid-state dc-ac 
inverters. Finally, because the operating generators’ frequencies are decoupled from the 
ac loads’ operating frequencies, the prime movers can function at their most efficient 
speeds. 

Some disadvantages to DC ZEDS must be acknowledged. The state-of-the-art 
solid-state switches used in power electronics today still limit the dc bus voltage to 
between 1500 volts and 2000 volts. Maximizing bus voltage is important in the effort to 
minimize conduction losses and therefore minimize conductor size. However, higher 
voltages pose concerns regarding system grounding and proper isolation of the dc-dc 
converters. An additional disadvantage is that designers need to more intelligently merge 
the various components within the system due to stability concerns. 

D. PREVIOUS DC ZEDS RESEARCH 

Several theses from the Naval Postgraduate School have investigated various 
aspects of DC ZEDS. These areas include modeling portions of the system to observe 
controllability and stability issues of both dc-dc converters and synchronous generators, 
as well as building and testing dc-dc converters and controllers [2-14]. Most recently, an 
investigation of frequency-based load sharing between current-mode-controlled buck 
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converters was conducted which included an rms frequency estimation algorithm and 
accompanying hardware to achieve the desired load sharing characteristics [15]. 

Previous work conducted in the fields of synchronous machine modeling and dc 
link filter stability issues were used as the foundation for this thesis. Modeling of an 
uncontrolled synchronous machine in a voltage-behind-reactance representation was 
conducted by Sudhoff and Wasynczuk [16] and Branson, Ciezki, Sudhoff, and 
Wasynczuk [17]. Research conducted by Sudhoff, Corzine, and Glover at the University 
of Missouri-Rolla, and Hegner and Robey at the Naval Surface Warfare Center makes 
observations on stability conditions in a simplified dc link model [18]. A technical report 
by Ciezki and Ashton at the Naval Postgraduate School discussed converter and converter 
control design methodology [19]. A dissertation by Belkhayat explored an impedance- 
based stability criteria for a stand-alone dc-dc converter [20]. This small-signal 
requirement was investigated by application of Nyquist stability criteria. Finally, 
research conducted by Belkhayat and Cooley at the Naval Surface Warfare Center, and 
Witulski at the University of Arizona investigated large-signal stability criteria for 
distributed systems with constant power loads [21]. 

E. THESIS GOALS 

This thesis explores various stability criteria for a representative ship-board dc 
power distribution system that includes a synchronous generator, a phase-controlled 
rectifier, a LC link filter, and a dc-dc converter. Specific limits on the filter parameters 
and converter input impedance which can be given as design specifications to contractors 
are the desired results. The analysis methods examined are a simplified constant-power 
dc link model, a model based on system impedances, and a Lyapunov-based, large-signal 
nonlinear representation. All models were studied using the Advanced Continuous 
Simulation Language (ACSL) together with specially written MATLAB script files. In 
addition, an ACSL model of the entire system was created for more advanced analysis. 
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F. 


CHAPTER OVERVIEW 


The derivation of modeling equations for a reduced-order synchronous 
machine/rectifier drive are detailed in Chapter II. The accompanying model code is then 
explained in detail. Chapter III examines a simplified dc link model and the limitations 
placed to ensure stability by the Hurwitz criterion. The results of simulations used to 
verify the limitations are presented. Chapter IV then examines stability criteria based on 
a system impedance model. Design of a dc converter is outlined, the required impedance 
values are determined and the Nyquist criterion is applied to the model. Nyquist plots for 
stable and unstable conditions are displayed. Chapter V covers the application of a 
Lyapunov-based large-signal stability method to a simplified system model to ensure that 
large system disturbances are evaluated in addition to the earlier small-signal models. 
Specifically, the equilibrium point criterion and mixed potential criterion are evaluated. 
Finally, Chapter VI presents the conclusions and discusses areas available for further 
research. 
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II. SYNCHRONOUS MACHINE - RECTIFIER MODELING 


A. REDUCED-ORDER MACHINE DEVELOPMENT 

The first step in developing an overall model is accurately modeling the electrical 
source, which in the system being evaluated is a three-phase, ac synchronous generator 
connected to a phase-controlled rectifier. A detailed, valve-by-valve simulation is 
computationally intensive and results in considerably long run times. Since the dynamic 
effects of the system are important for the analysis, a time-efficient reduced-order model 
is used for simulation and analytical purposes. Figure II-1 shows a synchronous machine 
depicted as a voltage-in / current-out model. The machine terminal voltages, therefore, 
are required as both inputs to the generator governing equations and for use in 
determining the commutation properties of the thyristors in the rectifier. Thus, there is an 
inherent algebraic incompatibility that must be resolved by the model synthesis. 

1. Approach 

The following methodical approach was used in the development of the reduced- 
order model. The machine stator electrical transients (known as the “fast” transients) 
were ignored for simulation purposes. A voltage-behind-reactance (VBR) representation 
of the synchronous machine was then derived. The VBR model was subsequently used 
.along with the different system modes to evolve an equation for the average direct current 
output voltage of the rectifier (V dc ). The source representation was completed by 
establishing the required synchronous machine stator currents. 
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2. Synchronous Machine Modeling Equations 


From reference [16], the equations used to describe a reduced-order synchronous 
machine are given. These equations were derived by setting the derivatives of the stator 
flux linkages per second equal to zero and fixing the reference frame speed equal to o) e 
(the electrical angular velocity). Further, assuming an isolated synchronous machine in 
this case, we find that o e = co r (the rotor electrical angular velocity). The stator voltage 
equations are as follows: 


\> r = 
V 


r r C + — 'FI 


s qs 


4 = r, 4 + 


co h 


(O h 


ds 


l}/r 


( 1 ) 

( 2 ) 


where v' s and are the machine terminal voltages expressed in the rotor reference 
frame (RRF), i r qs and are the stator currents in the RRF, ^ and are the stator 
flux linkages per second in the RRF, and co b is the base electrical angular velocity of the 
machine. The transformation matrix which relates the RRF quantities to the abc 
reference frame quantities is: 
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where r^, r u , and r fd are the referred winding resistances of the rotor coils. The 
referred values, designated by the primed quantities, are simply scaled versions of actual 
machine quantities. The scaling factor is related to the ratio of winding turns between the 


rotor and the stator (see [22].). 

Stator flux linkage per second equations are given by: 

+ <*) 

Rotor-referred flux linkage per second equations are expressed by: 

+ x ^t + 0 ( 9 ) 

~ X ikd i/d + X md (i^ +i u + if d ) ( 10 ) 

^fd = X,fd *fd + Xmd ifds + *kd + *fd ) (H) 
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where the various reactances are constants for a given machine provided operation 
remains magnetically linear. 

One expression for the torque may be written in terms of the stator currents and 
stator flux linkages per second: 

(i2) 

where P is the number of poles of the machine. 


The state variables for the system are x ¥' kd , m' fd , a r , and 0 r . The stator flux 
linkages per second are not state variables since we have ignored the derivative terms in 
those voltage equations. The next step is to express and in terms of state 

variables and stator current. Solving (10) and (11) for i M and i fd and solving (9) for , 
we find: 




-xj\ 


*kd 


1 

1 

p- 

'^-x m A 

}fd_ 

X'uX# ~X 2 md 



(13) 


% Y ' 
A kq 


X , 


X 


mq -r 
l qs 
kq 


(14) 


Equations (13) and (14) are substituted into the stator flux linkage per second 
Equations (7) and (8) to yield: 


as) 

where the q and d-axis subtransient reactances are denoted X q and X d , and the q and d- 
axis subtransient flux linkages per second are denoted 'Pj and %. It can be noted that 

the subtransient flux linkages per second are functions of the state variables of the 
synchronous machine. The next step is to then substitute Equations (15) and (16) into the 
stator voltage Equations (1) and (2) to give: 


10 



Y qs s qs 


v ds -'s'ds 


co. 


3 

CO 


It is customary to define: 


-xX+—% 

(17) 

, °>b 


b w b 

(18) 


(19) 



i 

k 

(20) 

<° b . 



then Equations (17) and (18) can be rewritten: 


®b 


v:,=r r c-^x:c+E: 

a b 


( 21 ) 

( 22 ) 


Equations (21) and (22) are the desired reduced-order voltage-behind-reactance (VBR) 
representation of the synchronous machine. 


3. Modeling Average DC Voltage 

The next step involved in modeling the electrical power source is to derive an 
expression for the average dc voltage using the VBR model together with boundary 
values associated with the conduction of the rectifier. In reference[16], the technique 
used was to find the average dc voltage during a single switching interval, since each of 
the six steady-state switching intervals in a three-phase rectifier are identical. During a 
single switching interval, the average dc voltage is a function of the phase voltages of two 
of the phases of die synchronous machine over the switching interval as shown by: 

e r =Y +/3 

y*-T7 f(Vfc-V.ye, (23) 

n/ i 
/ 3 &r=— + fi 
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where /? is the firing delay angle from the natural point of commutation with respect to 
the rotor position [16]. 

The VBR model can be used to show that: 


3 lO 8r^~~ + P 

n co h e r =-+p 


(24) 


The small effect of the stator resistance is ignored. 

The flux linkages per second can then be determined at the two limits using the 
DC currents at each limit to yield the following equation for the average dc voltage: 


V dc = 


3^3 


CO. 


3 o). 


n co b 


VCO 2 *^) 2 cos#- — 

n co b 


r K+ x 'S 


/*+sm^/!+ ^ 

V o 


(25) 


where a is the firing delay angle from the natural point of commutation with respect to 
the source electrical angle [16]. 


4. Calculation of Stator Currents 

The final step in modeling the power source is to calculate the stator currents. 
This is necessary since the referred rotor currents (Equations (13) and (14)) are expressed 
in terms of rotor-referred flux linkages per second and rotor reference frame stator 
currents. In the developed model, the rotor reference fra m e stator currents are needed to 
update the state of the machine during each time step iteration. The non-ideal case of the 
current flow through the six individual thyristors as a function of the three-phase 
synchronous generator output as shown in Figure II-2. As it would be expected, these 
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1 T 


/3 f*/> 

The non-ideal switching interval can further be divided into two distinct events. The first 
part of the switching interval, called the conduction period, consists of the section of time 
that there is only one device conducting current on each rail. The second part of the 
switching interval, called the commutation period, occurs when current is being 
transferred from one device to another device on the same rail. Thus: 

lqs,ave ^ qs.com ?qs,cond ( 28 ) 


i = i + i 

'ds.ave l ds,com l ds,cond 

If u is defined as the conduction angle, then Equations (26) and (27) can be 
subdivided into the conduction and commutation periods as follows: 


iqssom = —7 \i r qs (e r )de r 


3 §♦» 


3 P 

J<i (O r )dB r 


'3 f.j> 


1 T 

Ks,co nd =—? \iqs(0r)de r 


3 f+»+. 


j l P 

ids,cond = ~^T7 j Ids (fir @ r 


3 f +P +« 


Since the currents are known at the boundaries of the conduction interval, 
modeling equations for the conduction currents are easily obtained from Equations (32) 
and (33): 
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The only unknown in the above equations is the commutation angle (u). This can 
be found while examining the commutation period dynamics. Finding the commutation 
currents is a more difficult endeavor. During commutation, two of the phases are shorted 
together. As a result, the difference between the flux linkages per second of the two 
phases is a constant, specified in [16] as AT,. The quantity i as can be solved for in the 
abc reference frame yielding: 


K,-yf. 3T; cos|V +1 j + y; sin(e, + |j 

(A-;-x;)co S f2« r +|l+(x; + x;) 

V ' (36) 

(x; - X' d )cos(20 r + y) + \(K + X ',) u 

K - x' d )co/ w, +11+(x; + x;) 

where the aforementioned K x is: 

K,=S(- % cos p +>p; sin p)(x' d - X d )co/2/! + yj -i(x; + X \) /* 

(37) 

When commutation is complete, z as as a function of 0 r is zero. An implicit 
expression for the commutation angle u is found by setting i m equal to zero in Equation 
(36): 

/„.,(«)=*,-^[-^sinCS + ^+^cosCS + K)] 

- (x q -X d )cos\2p + 2u-^y±{X q +X d )l dc 
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Given the state of the machine, the dc link and the firing angle, the commutation 
angle can be determined by setting Equation (38) equal to zero and solving for u . This is 
done in the simulation by iteration to achieve a commutation angle within 0.001 radians. 
The iterative method used simply involves halving the range in which the commutation 
angle exists during each step until the range is narrowed to the above stated tolerance. 
Once u is calculated, the conduction current components can be determined using 
Equations (34) and (35). 

The q and d-stator currents in the RRF are defined during the commutation period 


as: 


Ks = 


4 = 


2V3 

3 

2V3 


r 


'as COS 


n 


^r+^J-4 sin ^r 


'as Sm 


e.+- 


+ I dc cos6 r 


(39) 


(40) 


V 

Equation (36) is substituted into the above equations, which are then substituted 
into Equations (30) and (31). The resulting integrals cannot be solved by analytical 
methods. Instead, a numerical method of integration called the Simpson method is used 
[16]. The Simpson method yields the following equations for the commutation currents 
over the commutation angle by utilizing a weighted measurement of the q and d-currents 
at increments through the commutation process, including the boundaries: 


qs,com 


L ds,com 


3 u 
xl2 


3 u 
n 12 






+4/1 


* ' ■'** Wf+f +2l ' gs +4, «* 


(41) 


4 ^ =>5+ - + ^4 L \ a _ a nu +4/* \ a _ e „ 3„ +/<* \ n o „ 

r 3 


0 r =j3 +—+- 
3 4 


3 2 3 2 


e r =p+-+u 


(42) 

Now that both the conduction and commutation currents are known, the average 
RRF stator currents are calculated with Equations (28) and (29). These currents are then 
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injected into the synchronous machine rotor flux linkage per second equations shown in 
Equations (15) and (16). 

B. DESCRIPTION OF ACSL CODE AND SIMULATION 

The derived model is programmed in the Advanced Continuous Simulation 
Language (ACSL). A listing of the .csl and .cmd files is documented in Appendix A. 
The .csl file contains the equation description of the synchronous generator-phase 
controlled rectifier in ACSL statement format. The .cmd file contains sets of runtime 
instructions to manipulate the model output into a usable form. ACSL was chosen for its 
execution speed and its strength and ease of working with and sorting continuous 
modeling equations. The feature is desirable due to the number of variables involved 
with the model as well as the interconnectivity of the numerous variables. 

1. .CSL File Description 

A block diagram representation of the model in terms of the required inputs and 
resultant outputs is shown in Figure II-3. The initial section of the .csl file holds the 
machine constants used in the test model. The assumed per unit electrical machine 
parameters were taken from reference [17] for the following: r s , X, s , X mq , X md , I dc , 

X' lkd , X\ fd , X lkq , r'w , r fd , r'^ and e xfd , where the leakage reactance terms with a 
subscript “7” are used to calculate the q and d-axis subtransient reactances (X q and X d ). 

The derivative section contains the equations that describe the system. Within the 
derivative section, there are two subroutine calls and two function calls. The two 
FORTRAN subroutines implement the Simpson method described in reference [16] to 
calculate the five different current components in equations (40) and (41). The first 
FORTRAN function calculates the commutation angle (u ) as described in an earlier 
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Figure II-3 Block Diagram of Synchronous Generator-Phase Controlled Rectifier 

ACSL Model 

section. The second FORTRAN function calculates the firing delay angle with respect to 
rotor angle (/?) using the same method as adopted for calculating u. 

The model accounts for the fact that /? must be set such that the thyristor is not 
gated until after it is forward biased by constraining /? to a minimum value, /? min , that 
denotes the earliest that a thyristor can be turned on. The firing angle for the thyristor 
with respect to the point its anode-to-cathode voltage becomes positive is then 
represented by a, which is the difference between p and fi min . 

2. .CMD File Description 

The .cmd file is configured to run a simulation to produce data points which 
define V dc and T e as a function of p . The resulting graphs for I dc = 0.9 pu are shown in 
Figure II-4. The effect of this is constant values for V dc and T e at values of fl less than 
the minimum value that represents gating the thyristor prior to it being forward biased. It 
can be seen from the graphs that fi min is at approximately 30°. With p equal to p mjn , 

V dc and T e are at their extreme per unit values of 1.4 and -0.85 respectively. V dc and T e 
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change linearly as J3 is increased in value when it is greater than /? min . V dc decreases 
from a positive value and T e increasing from a negative value, corresponding to generator 
action. When /? exceeds approximately 70°, the average dc voltage becomes negative 
and the average torque becomes positive, corresponding to motor action. This linear 
region lasts until /? reaches approximately 105°. 




Figure II-4 Beta vs Vdc and Te 

A detailed description of a model for the system source has now been presented. 
The next chapter discusses the evaluation of a simplified dc link-converter model for 
stability criteria in terms of major system parameters. 
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III. SIMPLIFIED DC LINK/CONVERTER MODEL 


A. DEVELOPMENT OF THE MODEL 

One characteristic predominant in dc-dc converters is the apparent small-signal 
negative input impedance which is a result of the high bandwidth regulation that makes 
the converter appear as a constant power sink. In Ref. [18], it is noted that one method 
that can be used to ensure stability is to increase the dc link capacitance. This is not a 
desirable remedy since it requires the use of more space and increases the weight 
contribution of the overall distribution system. In addition, since electrolytic capacitors 
have relatively low reliability, large banks are unattractive owing to the time-intensive 
task of identifying a shorted capacitor. The authors in [18] develop a simplified model, 
as shown in Figure III-1, to illustrate the cause of negative input impedance instability. 



Figure III-l Simplified DC Link Model 


Several simplifying assumptions are made. First, the prime mover and voltage 
regulator dynamics are assumed to have significantly long time constants. Second, all 
losses in the machine and rectifier axe assumed to be zero. Third, the rectifier is assumed 
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to be uncontrolled and finally, field voltage regulation maintains the source voltage equal 
to the commanded source voltage resulting in no need for subtransient information. 

The equations which describe the non-linear simplified model are then found to 
be: 


d K 

dt 


L 


= -—v dci -—i.+ 


K 

L. 


L Ves 


dv 


P' 


dci _ f j _ 

dt C e C e v dci 


(43) 


(44) 


where L e is a function of the filter inductance and the transient commutating inductance; 

R e is a function of the commutating inductance, machine stator resistance, and link filter 

resistance; and C e includes the filter capacitance and any input capacitance associated 

with the converter. P is the assumed constant power being drawn by the converter. 
Linearization of the above equations result in: 


d A i, 


1 


R 


- = —— Av rfd -—A i. + — Av. 


dt L : 
d Av. ; 


^ _ 1 A . P » 

7, 2 ^dci 

dt C e C e V dci 


(45) 

(46) 


Furthermore, the linearized converter input current is: 


Ai dci = 


dV.; 


P' 


dci L v dci J 


-P* 


Av dci = 2 


Av 


dci 


(47) 


v dci 


Thus, the small-signal converter input impedance may be found from: 


Z =^L- 


v dci 

tidci P' 


(48) 


Therefore, for a given power being supplied by the converter and a given input voltage 
v dcj , the input impedance must be negative. It is this negative resistance which sets up 
the possibility for oscillations with the input LC filter. 
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B. EVALUATION OF THE MODEL 

Confirmation that the negative input impedance results in system instability is 
obtained by examining the characteristic equation of the model. In order for the system 
to be stable, a necessary and sufficient condition of the Hurwitz criterion requires that all 
coefficients of a second-order characteristic equation be positive. From Equations (45) 
and (46), the characteristic equation for the system is: 


4 C e v 2 dci 


44 C e L e v dd 


The resulting stability requirements are then: 


R e P* 

>-r 

L. C.vi 


yr>R e ( 51 ) 

Equation (51) is typically satisfied by matching the load characteristics to the given 
source resistance. Thus, Equation (50) becomes a first-order design criterion for the link 
filter parameters. A rearranged version of Equation (50) illustrates that we can ensure 
small-signal stability by including additional bus capacitance such that: 

C P* 

— >-V (52) 

4 RA* 

ACSL code for the .CSL and .CMD files which model both the linearized and 
non-linearized representations of the model is given in Appendix C. Simulations were 
conducted in which P* was assumed to be 200 kW, v dcj was assumed to be 1100V, R e 
was assumed to be 0.01 Q , and L e was assumed to be 0.5 mH. This resulted in requiring 
C e to be 8265 ju F to just be stable in the event of large disturbances. 

The simulated output of a transient from 10% rated power to 100% rated power at 
t im e 0.5 seconds and a return to 10% rated power at time 3 seconds is shown in Figure 
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III-2. A value of C e = 8500 juF was used to better show the stable response due to the 
low damping of this configuration. The system can be seen to be stable after both the up- 
power and down-power transients. The simulated output of a similar transient sequence 
with the down-power transient occuring at time 0.6 seconds is shown in Figure III-3. For 
this simulation, C e was lowered to 826.5 jjF . It is shown that following the up-power 

transient, the system becomes unstable. Stability is not recovered following the down- 
power transient. 

We have now established the governing stability criteria dictated by a simplified 
dc lmk/converter model. In the next chapter, we examine a more detailed model that 
defines the requirements for stability using subsystem impedances. 



Figure III-2 Up-Power, Down-Power Transient, C £ = 8500 MicroFarads 
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IV. IMPEDANCE-BASED STABILITY CONSTRAINTS 


A. STABILITY CRITERIA USING SUBSYSTEM IMPEDANCES 

One analysis that can be conducted to determine the stability of the overall system 
involves determining the input impedance of the dc converter and the output impedance 
of the filter. Together, the poles of the overall system transfer function are found. The 
transfer function can then be examined to determine if the Nyquist stability criteria are 
met. The steps involved in this process are: design the converter to be used, determine 
the converter controller gains, calculate the converter input impedance, calculate the filter 
output impedance, and use the results to analyze stability. For analysis purpose, two 
converters are developed for use in simulations. The first is designed for a low switching 
frequency used for a high-power hard switching process. The second will be designed for 
a higher switching frequency used for a high-power soft switching process. 

B. DESIGN OF THE DC CONVERTER 

The layout of a common step-down dc-dc converter in a Buck Chopper 
configuration is shown in Figure IV-1. To design the converter, the first issue which 
must be resolved is the system requirements. To establish a baseline for the design, we 
require maintaining continuous current conduction mode for all loads that are greater than 
ten percent of the rated power. For this thesis, the design specifications listed in Table 
IV-1 are assumed. 
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Specification 

Value 

Rated Power (P rat ) 

200 kW 

Filter Output Voltage (V dc ) 

1100 V 

Load Voltage (V out ) 

950 V 

Ripple Voltage (A V out ) 

5 V 

Hard Switching Frequency (f shard ) 

4 kHz 

Soft Switching Frequency (f s soft ) 

20 kHz 


Table IV-1 Design Specifications for the DC-DC Converter 



R (load) 


Figure IV-1 DC-DC Converter - Buck Chopper 
1. Design of the Components 


The authors in [19] outline a sound methodology for converter design. The 
nominal duty cycle (D 0 ) of the switch is defined as the percentage of time that the switch 
is closed. For the designed converter, assuming continuous current conduction, we get: 


Z> = 



950F 

1100F 


0.8636 


(53) 
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Next, based on the rated power and the desired output voltage, the rated load resistance 
can be found by: 


R 


rat 


VL foor ) 2 

P rm 200kW 


4.5125 Q 


(54) 


The critical resistance (R cri ,) is the load resistance at the boundary between continuous 
and discontinuous current conduction modes. For this converter: 

R cril =10 xR ral = (l0)x (4.5125Q)= 45.125Q (55) 

The main components of the converter can now be selected. The hard switching 
converter parameters are determined first. The critical inductance (l c „,) is the minimum 
inductance required to maintain continuous current conduction mode during all 
conditions required by the design specifications. The critical inductance is found by: 

4*,=ff(>-A,)=7 69AmH (56) 

s 

In order to maintain a comfortable margin, L is assumed to be ImH. 

The minimum capacitance is designed around ensuring that the steady-state 
output voltage ripple specification is met: 

C„i.= W fc-0 = 2°2- 4 /<r (57) 

Once again, a margin is included to ensure that the capacitor is large enough to provide 
sufficient energy during a transient to maintain output voltage. The capacitor selected 
should not be too large, however, because this would result in small changes in output 
voltage causing the duty cycle to oscillate between full-on and full-off. C is assumed to 
be 2mF. 


2. Converter Controller Design 

The buck chopper exhibits an oscillatory response as a result of the undamped 
L-C filter at the converter output. To ensure that the converter operates satisfactorily 
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during transients and to provide steady-state converter regulation, a multi-loop 
feedforward controller is utilized for duty cycle control. The governing controller 
equation is [19]: 

4 ) = V -f- +- v„ )+ h„K, j(v r , f - v„ )dt + h,K,(h - i m ) 

V dc 


(58) 

where 1^, h„, and hj are the feedback gains, K,, and Kj are the voltage and current sensor 
gains, V ref is the desired output voltage, and i out is the sensed output current. This 
controller design assumes that the converter is operating autonomously. Additional terms 
are required to parallel multiple dc-dc converters. The first component of Equation (58) 
provides feedforward action to ensure rapid response to disturbances in the input voltage. 
The second and third terms implement a proportional-plus-integral (PI) outer voltage 
control while the fourth term realizes an inner current control. Effectively, this last 
component represents the feedback of the derivative of the capacitor voltage so that a 
PID-type contol is realized. 

By assuming that the desired output voltage (V K/ ) remains constant and the 

converter load is purely resistive in nature, we find that the operation of the converter can 
be described by a set of averaged linear ordinary differential equations with the following 
characteristic polynomial: 


s 3 +d 2 s 2 +d,s + d 0 (59) 

where d 2 , d l5 and d 0 (functions of the converter parameters and control gains) describe the 
system pole locations and, ultimately, the converter’s transient response. Upon selecting 
desired characteristic polynomial coefficients, we establish gain relationships as follows: 


KK,=y- 

V dc 


4 - 


LC 


hK = 


d n LC 


dc 


(60) 

(61) 
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(62) 


« = 


-L 


V 


dc 


<L~ 


RC, 


The utilization of a third-order Bessel function polynomial as a system prototype 
is known to result in a desirable system response. In addition, it minimizes the distance 
that the open-loop poles must be moved and hence minimizes the control effort. The 
specific poles for the prototype system are located at: 

5 , 2 = -0.7455 a>*/_ jOJl 12 co (63) 


s 3 = -0.9420 co 


(64) 


where co is the desired closed-loop bandwidth. The radian frequency {co) must meet 
three criteria to provide optimized results. The first criteria is that co must be at least one 
decade below the switching frequency to prevent undesirable controller-switch 
interactions. The second criteria is that co should be minimized to prevent additional 
noise from being introduced into the system through the controller’s operation. The final 
criteria is that co must be large enough to achieve the desired settling time following a 
transient. As a result of these criteria, co was chosen for the hard switching converter to 
be 1500 rad / sec . The poles for the designed converter are then: 

s 12 =-1118.25 Y_j 1066.8 


= -1413 


The system characteristic equation is then: 

5 3 + 3639.5s 2 +5526354.85 + 3351129059.4 
and the gain selections become: 


h v K v = 0.00914 


h n K v = 6.093 


(65) 

( 66 ) 

(67) 

( 68 ) 

(69) 

(70) 


h i K i = -0.003299 

Similar calculations were made for a soft switching buck chopper with a 
switching frequency of 20 kHz. L crU was 153.9 fj. H, with a value of 200 /u H chosen to 
provide a margin. C min was found to be 40.5 n F, with a value of 400 jj. F chosen to 
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ensure adequate energy storage capacity. The desired system poles utilizing the Bessel 
prototype with co = 3000 rad / sec are: 


j u =-2236.5 y_j2 133.6 (71) 

= -2826 (72) 

The resultant gains are: 

h v K v = 0.00070508 (73) 

h n K v = 1.9636 (74) 

h,K, =-0.001317 (75) 


3. Filter - Converter Model in ACSL 

The .csl and .cmd files for an ACSL model of a dc link filter attached to the input 
of a dc-dc converter with feedback control are contained in Appendix D. The model 
specifically contains the component parameters derived for the hard switching converter. 
It assumes an ideal dc input voltage. The initial conditions for the simulation are that the 
filter capacitor is fully charged, and the converter is secured. Converter operation is 
initiated by ramping up the reference voltage. Converter control is implemented with the 
method outlined in the previous section. 

Several simulations were run with various filter parameters. The behavior of the 
system with changes in the filter capacitance were of particular interest. When filter 
parameters were within those specifed by the stability requirements, the output was stable 
as expected. However, of interest was the observation of a region of stability outside of 
the criteria bounds established in the stability analysis. The band was observed to exist 
when filter capacitor values between approximately 25 /iF and 110 juF were used in the 
simulations. It is believed that this region of stability is a result of the feedback control 
for this particular set of parameters. In general, changes in the filter capacitance resulted 
in changes in the ripple voltage magnitude and frequency. 
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C. IMPEDANCE CALCULATIONS 

A representation of a source-filter-converter-load is shown in Figure IV-2, as also 
shown in reference [20]. The overall transfer function of this system is shown to be: 

v, l + Z./Z, 

where Z s is the output impedance of a line LC filter and source, and Z, is the input 
impedance of a loaded converter with closed-loop control. Stability of the system can be 
determined by analyzing the poles of this transfer function. This section will describe 
how Z s and Z, are determined for this system as well as analyze the behavior of the 
poles of the system. 
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Figure IV-2 DC-DC Converter With Input Filter 
1. Determination of Converter Input Impedance 


Due to the non-linear nature of a buck chopper, determination of the input 
impedance needs to be performed on a linearized version of the model. Referring back to 
Figure IV-1, we find the input current is related to the averaged inductor current by: 

h =di L (77) 

where^ is the converter duty cycle. Linearization of Equation 67 yields: 

Ai b =d 0 bd L + i u Ad (78) 


33 





where Ad is the perturbation variable given by : 

Ad = d-d 0 (79) 

The constants d a and z io are operating point values for those variables. By assuming 
AV ref is zero, the linearization of Equation 48 can be shown to yield: 

Ad = K] AVj n + T x {s)Av c (80) 

where AT, is: 


K: = 


- V: 


ref,o 


(81) 


m,o 


and T x (s) is: 


T ( 5 ) = s fyi 


(82) 


Further definitions that are made involve the transfer function between the variation in 
input voltage (Av m ) and the variation in the state variable inductor current (A i L ) as well 
as the variation in the state variable capacitor voltage (Av c ). These relationships are: 

Ai L =T 2 (s)Av in (83) 

Av c =T 3 (s)Av m (84) 

The MATLAB code used to determine the coefficients for T 2 (s) and T 3 ( 5 ) is 
contained in Appendix E. The method used to find these involves deriving the state- 
space form of the linearized system equations, then utilizing the state-space to transfer 
function (ss2tf) MATLAB function with a C matrix that zeros all outputs except A i L for 
the T 2 (s) calculation and Av c for the T 3 (s) calculation. 

The desired input impedance is defined by: 

Av,, 


Z f = 


Ai k 


(85) 


Equations (68), (70), (73), and (74) are then substituted into (75) to obtain: 
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( 86 ) 


JMs) + li,K,+i u T,(s) T,(s) 

The coefficients for the s terms are also determined by the MATLAB code in Appendix 
E. Input impedance calculations are performed for both the hard switching and soft 
switching converter options. Note, Equation (86) holds for the specific multi-loop 
control with feedforward under consideration. 

2. Determination of Filter Output Impedance 

Figure IV-3(a) is a representation of the simplified network looking back through 
the filter to the source. The network which is used to find the Thevenin equivalent 
impedance is shown in Figure IV-3(b). By appropriately combining series and parallel 
impedances: 



(a) (b) 

Figure TV-3 Filter Output Impedance (a) Simplified System Diagram, (b) Thevenin 

Equivalent Impedance Network 
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3. Impedance-Based Stability Criteria 

The type of stability used for this particular analysis is Bounded Input Bounded 
Output (BIBO) stability. Per reference [23], “A system is stable if every bounded input 
excites a bounded output. Otherwise the system is said to be unstable.” A system is 
determined to be stable if no poles of the open-loop transfer function are located in the 
open right-hand s-plane. Referring back to Equation 76, it is known that h c is designed 
to be stable through the placement of poles far in the left-hand s-plane with the effective 
PID (multi-loop) controller. The filter transfer function (h f ) is also always BIBO stable 

(it is a strictly passive network); therefore, the product of these two terms will conversely 
be BIBO stable. As a result of this, stability of the overall system can be determined by 
analyzing the denominator of Equation 76 using the Nyquist stability theorem. In order 
for the system to be stable, the Nyquist stability theorem requires that the Nyquist plot of 
the system encircle the (-1 + Oj) point counter-clockwise the number of times equal to the 
number of poles in the open right-hand s-plane. Since there are no system poles in the 
open right-hand s-plane in this case, no encirclements of the (-1 +0j) point should occur. 
To further ensure that stability exists, the Nyquist plot of the system is confined to the 
unit circle: 

<1 for - oo < &) < co (88) 

With these requirements in place, BIBO stability is guaranteed. 

4. Evaluation of the Analyzed System 

A MATLAB file was written to examine the stability of the system using various 
filter parameters. This code is included in Appendix E. The Nyquist plot using the hard 
switching model together with filter parameters L flUe = 0.5 mH and C fiUer = 8265 fiF, 

and a source resistance estimated at R source = 0.01 Q is shown in Figure IV-4. It can be 



36 



seen that the Nyquist contour is positioned within the unit circle. An examination of the 
poles shows that there are no open right-hand s-plane poles. It is then confirmed to be 
BIBO stable. 

Nyquf»1 Diftoram* 



Figure IV-4 Nyquist Plot of BIBO Stable System 


The result of lowering the filter capacitor to 6000 juF is shown in Figure IV-5. 
The Nyquist contour encircles (-1 +0j) twice now. There are still no right-hand s-plane 
poles; therefore, BIBO stability is not assured. 

A noteworthy anomoly was encountered as a result of the method used to 
calculate the controlled converter’s poles and zeros. There exists a pole-zero cancellation 
that initially seemed not to occur. In fact, it appeared that the system had a RHP pole. 
Further investigation revealed that the cause was roundoff errors induced over multiple 
intermediate steps in the calculations. 

To this point, two small-signal models have been evaluated to establish stability 
criteria for the dc link filter. In the next chapter, Lyapunov-based large-signal stability 
contraints are applied to a basic system model. 
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Nyqulst Diagrams 



Figure IV-5 Nyquist Plot of Non-BIBO Stable System 
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V. LYAPUNOV-BASED LARGE-SIGNAL STABILITY CONSTRAINTS 


The previous stability criterion, specifically the small-signal models, apply to 
small disturbances in the system. In Reference [21], the authors examine a stricter set of 
criteria which act to bound the source resistance for a given set of filter, source voltage, 
and load power parameters. They note that when small-signal models are used, the 
stability margins that are used for the final design are forced to be very conservative to 
ensure stability during a large disturbance. The two criteria specifically examined are the 
equilibrium point criterion and the mixed potential criterion. This thesis then applies 
these criteria to representative values being considered for the DC ZEDS project. 


A. EQUILIBRIUM POINT CRITERION 


The system examined is the same as in Figure III-1. Equations (43) and (44) are 
the equations that describe the states of the system. The equilibrium points of these 
equations are found by setting the derivatives of i r and v dcj with respect to time to zero: 


and 



(89) 



(90) 


where capital letters are used to distinguish steady-state quantities from instantaneous 
values. 


In order for the above equilibrium points to be stable, they must be positive and 

(vy 


real. For the v dci operating point to be real, R e P* must be smaller than the 
This results in the condition: 


term. 


v ^ y 
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V 2 

es 


( 91 ) 


R e <—r 
e 4 P 


The authors in [21] argue that since this is a hard requirement, it is prudent to apply a 
safety factor of two when using this criteria. The combined source-filter resistance must 
satisfy: 


R 


e,m arg 



(92) 


B. MIXED POTENTIAL CRITERION (SINGLE CONSTANT POWER 
LOAD) 


The authors next examine the mixed potential criterion, which was originally 
developed specifically for non-linear systems. The following mixed potential equation is 
found using the same system as that used for the equilibrium point criterion: 

u = -\ R e' 2 r + \— dv dci + *V( v « ~ v dd) (93) 

Z V dci 

Lyapunov-based stability conditions are then placed on this system which results in the 
following requirement to ensure that the system response to a large disturbance is not 
oscillatory [see 21]: 



Since this is not a hard value, no safety margin is required for this criteria. Thus, the 
design specification on the source plus filter resistance combines Equations (92) and (94) 
to yield: 



(95) 


Thus, a low equivalent source resistance is attractive from the standpoint of 
regulation and efficiency, but it in turn may require a large C e to ensure system stability. 
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The limitation of this analysis is that a fairly simple representation of the converter is still 
employed. In particular, little can be deduced as to the impact of control bandwidth on 
the stability conditions. Also, it relies on a very basic representation of the synchronous 
machine and rectifier. Since the criteria does not guarantee a stability margin, simulation 
must be used to assess transient response characteristics. The attendent simplicity, 
however, does provide a convenient, simple tool which can be used to provide baseline 
conclusions. 


C. MIXED POTENTIAL CRITERION (MULTIPLE CONSTANT POWER 
LOADS) 


In the dc power distribution system envisioned for DD-21, there are multiple 
converters and loads, each with their own link filters, operating in parallel supplied by the 
same source. This may represent a generator-rectifier supplying a dc bus which, in turn, 
delivers power to the various zones within the ship through separate regulated dc-dc 
converters with associated input filters. It is important to understand how the additional 
converters affect the constraints placed on the filter design. The system representation of 
a multiple constant power load system is shown in Figure V-l. The operating point 
criterion is applied to this configuration to yield: 


V 2 

^_ _es _ 

* 4(P\ + P2 + P3) 


(96) 


For simplicity, the case is considered where each of the filters and constant power 
loads are identical, therefore L x =L 2 =L 3 = L , C, = C 2 = C 3 = C , and 
P\ = P2 = P3 = P. Then the Lyapunov stability criterion requires establishing bounds 
on the norm of [see 21]: 


K = 


-4l 


R e Jc 


-2 2 

1 -2 
0 1 


(97) 


41 



If an upperbound of the norm of K is constrained to be less than unity for 
stability, the following relationship for source resistance holds: 

R e « 2 -n 

Vc 

where n is the number of converters. For the case of three converters illustrated in 
Figure V-l: 

is 



Figure V-l Multiple Constant Power Load System Representation 



( 100 ) 



The ACSL code of a single load system is contained in Appendix F. Several 
simulations of the system were conducted to verify the validity of the mixed potential and 
equilibrium point criteria. The values used for the simulation were P' = 200 kW, v es = 

1100 V, L e = 0.5 mH, and C e — 8265 fi F. These parameters result in boundaries for 

R between 0.2460 Q and 1.5125 Q. When an R of 0.017 Q is used, which is below the 
lower limit defined by the mixed potential criterion, it can be seen in Figure V-2 and 
Figure V-3 that the system is stable for a small 5% disturbance in the system voltage but 
is unstable for a large 25% disturbance in the system voltage. In both figures, the 
disturbance is introduced at time 0.5 seconds. The initial transients are due to the system 
reaching steady-state from the initial simulation conditions, which were approximated to 
be close to the anticipated steady-state values. 



T lime* 


Figure V-2 System Response With R=0.017 Ohms, 5% Drop 
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Figure V-3 System Response With R=0.017 Ohms, 25% Drop 


Likewise, when a value for R of 0.26 Q is used, it can be seen in Figure V-4 and 
Figure V-5 that the system is stable for both small (5%) and large (25%) disturbances in 
the system voltage. Like the previous simulations, in both figures, the disturbance is 
introduced at time 0.5 seconds. The initial transients are again due to the system reaching 
steady-state from the initial simulation conditions, which were approximated to be close 
to the anticipated steady-state values. 



Figure V-4 System Response With R=0.26 Ohms, 5% Drop 


44 




Load I Current Load 1 Voltage 





46 



VL INTEGRATED SYSTEM REPRESENTATION 


A. USEFULNESS OF AN INTEGRATED SYSTEM 

The models used to analyze the stability criteria were all limited in scope to 
varying degrees due to necessary simplifying assumptions. A more detailed model of the 
system under consideration would be useful to verify the validity of stability criteria 
established by different methods. Models of both the synchronous generator-phase 
controlled rectifier pair and the filter-converter pair have been presented earlier in this 
thesis. While useful, these models are also limited in that the dynamics of the load are 
not experienced in the generator-rectifier model and the dynamics of the source are not 
experienced in the filter-converter model. 

B. INTEGRATED SYSTEM MODEL 

A circuit diagram of the system modeled by the ACSL file contained in Appendix 
G is shown in Figure VI-1. A block diagram of the model which clearly shows the 
various subsystems represented in the model and the interaction of the system parameters 
between the subsystems is shown in Figure VI-2. The four main subsystems represented 
in the model are the synchronous machine and rectifer, the generator field regulator, the 
link filter and the dc-dc converter. The synchronous generator/rectifier section is based 
on the equations derived in Chapter II. With a dynamic system model, the need for 
voltage control of the generator becomes evident. The derivations for the equations 
which govern the synchronous machine voltage regulation are contained in Appendix B. 
The link filter and dc-dc converter section is identical to that discussed in Chapter IV. Of 
interest is the fact that the synchronous generator/rectifier/regulator portions are modeled 
in per unit parameters, while the filter/converter portions are modeled in actual 
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component parameters. This is accounted for in the code by utilizing conversion factors 
on the system variables passed between the two sections. 


Idc 4 



Figure VI-1 Circuit Diagram of Integrated System 



Figure VI-2 Block Diagram of Integrated System 







Several simulations were run using the model to verify the validity of both the 
model and the design criteria evaluated in this thesis. The results of an up-power 
transient from 10% rated power to 100% rated power at time 0.52 seconds and the 
subsequent down-power transient back to 10% rated power at time 0.57 seconds are 
shown in Figure VI-3. The system’s component values were all such that the various 
stability requirements were met. It is shown that the system response is stable and has an 
acceptable response time. 



Figure VI-3 Up-Power, Down-Power Transient, Cf = 8265 uF 
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The same up-power, down-power transient with the filter capacitor value changed 
to 400 juF is shown in Figure VI-4. While this value of capacitance is enough to 
maintain stability at the 10% rated power level, it is not sufficient to maintain stability at 
100% rated power. This is confirmed by the stability criterion presented in earlier 
chapters. Large, sustained oscillations are obtained on the dc bus. It is important to note 
that the dc-dc converter controller maintains the load voltage at the desired reference 
voltage even with the sustained oscillations present at the converter input. 



QO 



Figure VI-4 Up-Power, Down-Power Transient, Cf = 400 uF 
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A couple of nuances that were encountered while conducting these system studies 
deserve further investigation. First, optimization of the numerical integration time step 
should be performed to uncover at what point the system becomes numerically unstable 
and at what point it becomes inaccurate. Further, the steady-state duty cycle did not 
appear to settle at the anticipated value of 0.88. It seemed to hone in on a value 
approximately 5% higher. This anomoly should be explained. Third, the filter current 
should be bounded so that there are no occasions when it can go negative. Finally, the 
voltage/current scaling procedure should be replaced by a more methodical scheme. 

The ACSL model of the representative system has been shown to be a useful tool 
to verify the results of stability analyses conducted by any method. The next chapter will 
present a s ummar y of the findings of this thesis and discuss future work that can be 
performed to futher knowledge about this issue. 
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VII. CONCLUSIONS 


A. SUMMARY OF FINDINGS 

Two small-signal models and a large-signal model of a representative dc electric 
power distribution system are examined to determine the design criteria for the dc link 
filter in the system. The Lyapunov-based large-signal stability criterion is found to be the 
most conservative of the three analyzed. The next most conservative is the simplified dc 
link model with the Hurwitz criteria for stability applied. The least conservative is the 
application of the Nyquist criteria to the impedance-based model. Simulations of the 
different simplified models were used to illustrate the validity of the filter design 
limitations defined by each criterion. 

When component values that meet the simplified model stability criteria 
examined are used, detailed model simulations show that these criteria are conservative 
enough to be considered valid for hardward design. A supplemental finding worthy of 
comment is the region of stability outside of the criteria bounds established upon 
investigation of the detailed filter-converter model. The region exists as a result of the 
feedback controller utilized by the converter to control the duty cycle. 

B. FUTURE WORK 

This thesis motivates the investigation of other more detailed aspects of stability 
concerns in a DC ZEDS design. First, further investigation into the stability concerns 
involved with operating multiple controlled converters in parallel can be undertaken. The 
simplified models used in this thesis ignored the effects of controlled converters. The 
relationship between control bandwidth and stability should be investigated. An 
additional issue is raised when the load on the dc system is an inverter vice a converter. 


53 




The input characteristics of a tightly-regulated inverter differ from those of a dc-dc 
converter. 

This thesis has focussed on meeting the system stability requirements by proper 
selection of the passive components in the device input filters. Further study should be 
performed with regards to introducing nonlinear compensators into the system. These 
compensators would preclude the need for large capacitor banks which are both difficult 
to troubleshoot and have marginal reliability. Finally, the active compensation 
techniques may be applied to the synchronous machine field regulator, the bus voltage 
phase-controlled rectifier, or with regards to the dc : dc converter controls. Many research 
opportunities are currently available in this area and should be explored. 

Factors that were assumed negligible in this thesis that may need to be examined 
in the future include the effects of line inductance and parasitic capacitance in a full size 
system. Additionally, the effects of electro-magnetic interference (EMI) due to the high 
switching frequencies used in the converters was not considered. 

In conclusion, various stability criteria have been examined to address the 
stability issues involved with having link filters as an integral part of a dc electrical 
distribution system. The analysis and tools crafted in this thesis form a strong foundation 
for establishing design specifications for the L-C filters in the DC ZEDS program. 
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APPENDIX A: SYNCHRONOUS MACHINE - RECTIFIER ACSL MODEL 


A. .CSL FILE 


PROGRAM !"tues.csl" 


1"Greg Greseth, LT, USN - 8 Jan 99" 

! "Last Revised - 15 Feb 99 

•"This ACSL file simulates a reduced order model for a 


synchronous" 
are" 


! "machine and three phase rectifier. 
!"defined in 'per unit' values." 


The machine parameters 


INITIAL !"Section in which constants and program 


parameters 

are defined" 



ALGORITHM ialg = 5 

!"Runge-Kutta 4th order 


CONSTANT tstop =0.5 

!"stop time for 

integration 

r? 



MAXTERVAL maxt = l.e-5 

!"max integration step 

size" 


MINTERVAL mint = l.e-6 

!"min integration step 

size" 


CINTERVAL cint = 5.e-4 

!"data communication 

interval" 

CONSTANT wb = 377. 

CONSTANT rkgpri = 0.0613 
CONSTANT rkdpri = 0.24 

CONSTANT rfdpri = 0.00111 
CONSTANT J = 1. 

CONSTANT T1 = 1. 

CONSTANT B = 0. 



CONSTANT P = 2. 

CONSTANT Xls =0.08 

CONSTANT Xmq = 1.00 

CONSTANT Xlkqpri = 0.330 

CONSTANT Xmd = 1.77 

CONSTANT Xlfdpri = 0.137 

CONSTANT Xlkdpri = 0.334 

CONSTANT rs = 0.00515 

!"number of poles" 


CONSTANT exfd = 1. 

!"CONSTANT beta = 0.523599" 

!"CONSTANT alpha =0." 

CONSTANT alfahat =0.0 

!"2.15" 


CONSTANT del = 0.0 

CONSTANT PI = 3.141593 

CONSTANT radl20 = 2.094395 
CONSTANT rad60 = 1.047198 
CONSTANT rad30 = 0.523599 

! " 0.9 " 
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u = 0.05 !"initialize u for use in first findu function 

call" 

’"define the subtransient reactances" 

Xqpripri = Xls +(Xmq*Xlkqpri)/(Xlkqpri+Xmq) 

Xdpripri = 

Xls+(Xmd*Xlfdpri*Xlkdpri)/(Xmd*Xlfdpri+Xmd*Xlkdpri& 

+Xlkdpri*Xlfdpri) 

Vfdpri = (exfd*rfdpri)/Xmd 

COMMON / forsub / 

rad30,rad60,radl20,Xqpripri,Xdpripri,del,PI 
END !"of initial section" 

DYNAMIC 

TERMT(t.ge.tstop, ’Termination on time limit 1 ) 

DERIVATIVE 

!"state equations for psikqpri, psikdpri, psifdpri, wr 

and thetar" 

psikqpri = INTEG(-wb*rkqpri*ikqpri, 0.0) !"0.049" 

psikdpri = INTEG(-wb*rkdpri*ikdpri, 1.0) !"0.293" 

psifdpri = INTEG(-wb*rfdpri*ifdpri+wb*Vfdpri, 1.077) 

!"0.460" 

!"wrm = INTEG(((1./J)*(Te-Tl-B*wrm)) , 0.)" 

!"relationship between wrm and state wr 

!"wr = (P/2.)*wrm" 
wr = 377. 

thetar = INTEG(wr, 0.) 

’"thetar = mod(theta,rad60) + rad60d" 

!"electrical torque equation" 

i »»Te = ( (3. /P) / (4 . *wb) ) * (psidsr*iqsr-psiqsr*idsr) " 

Te = psidsr*iqsr-psiqsr*idsr 

’"solve first for ikdpri, ifdpri and ikqpri" 

ikdpri = (1. / ( (Xlkdpri+Xmd)*(Xlfdpri+Xmd)- 
Xmd**2.))*((Xlfdpri+Xmd)& 

*(psikdpri-Xmd*idsr)-Xmd*(psifdpri-Xmd*idsr)) 
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ifdpri - (l./((Xlkdpri+Xmd)*(Xlfdpri+Xmd)-Xmd**2.))*( 

Xmd*(psikdpri& 

-Xmd*idsr)+(Xlkdpri+Xmd)*(psifdpri-Xmd*idsr)) 

ikqpri = psikqpri/(Xlkqpri+Xmq)- 
(Xmq*iqsr)/(Xlkqpri+Xmq) 

!"equations to eliminate psiqsr and psidsr" 

psiqsr = Xqpripri*iqsr+psiqpripri 

psidsr = Xdpripri*idsr+psidpripri 

!"define the subtransient flux linkages" 

psiqpripri = (Xmq/(Xmq+Xlkqpri))*psikqpri 

psidpripri = 

((Xmd*psikdpri)/Xlkdpri+(Xmd*psifdpri)/Xlfdpri)& 

/(1.+Xmd/Xlfdpri+Xmd/Xlkdpri) 

'"define Eqpripri and Edpripri" 

Eqpripri = (wr*psidpripri)/wb 

Edpripri = -(wr*psiqpripri)/wb 

{"define phi to use to determine alpha" 

phi = ATAN2(Edpripri,Eqpripri) 

!"call function to determine beta" 

beta = findbeta(psidpripri,psiqpripri,alfahat) 

!"define alpha" 

alpha - beta - phi 

!"desired reduced order voltage-behind-reactance 

equations" 

Vqsr = rs*iqsr+(wr/wb)*Xdpripri*idsr+Eqpripri 

Vdsr = rs*idsr-(wr/wb)*Xqpripri*iqsr+Edpripri 

!"equation for average DC voltage out of rectifier" 

Vdc = (3.*sqrt(3.)/PI)*(wr/wb)*(sqrt(psiqpripri**2& 
+psidpripri**2))*(cos(alpha))- 
(3./PI)*(wr/wb)*(((Xqpripri+Xdpripri)/2.)*dcl& 

+(Xdpripri-Xqpripri)*dcl*sin(2.*beta+rad30)) 

{"need to define K1 for use in subroutines" 

Kone = sqrt(3.)*{- 

psiqpripri*sin(beta)+psidpripri*cos(beta))& 
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+((Xdpripri-Xqpripri)*cos(2.*beta+radl20)-0.5*(Xqpripri+Xdpripri))*dcl 

!"call function to determine u using ’half“it 1 method" 
u = findu(Kone,psiqpripri,psidpripri) 

!"now with the new u, we can determine conduction 

currents" 


-cos(beta+u+rad60 

-sin(beta+u+rad60 

method" 


iqsrcond = ((2.*sqrt(3.))/PI)*dcl*(cos(beta+radl20)& 

) ) 

idsrcond = ((2.*sqrt(3.))/PI)*dcl*(sin(beta+radl20)& 

)) 

!"calculate the commutating currents using Simpson 
!"first, define the 5 different values of thetar to be 


used" 


thetarl = beta+rad60 
thetar2 = beta+rad60+u/4. 
thetar3 = beta+rad60+u/2. 
thetar4 = beta+rad60+(3.*u)/4. 
thetarS = beta+rad60+u 

["call each q and d subroutine 5 times to calculate 

components" 

PROCEDURAL (Q1 = thetarl,Kone,psiqpripri,psidpripri) 
CALL 

simpsonq(thetarl,Kone,psiqpripri,psidpripri,Q1) 

END !"of procedural" 

PROCEDURAL (D1 = thetarl,Kone,psiqpripri,psidpripri) 
CALL 

simpsond(thetarl,Kone,psiqpripri,psidpripri, D1) 

END ["of procedural" 

PROCEDURAL (Q2 = thetar2,Kone,psiqpripri,psidpripri) 
CALL 

simpsonq(thetar2,Kone,psiqpripri,psidpripri, Q2) 

END ["of procedural" 

PROCEDURAL (D2 = thetar2,Kone,psiqpripri,psidpripri) 
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CALL 

simpsond(thetar2,Kone,psiqpripri,psidpripri, D2) 

END !"of procedural” 

PROCEDURAL (Q3 = thetar3,Kone,psiqpripri,psidpripri) 
CALL 

simpsonq(thetar3,Kone,psiqpripri,psidpripri,Q3) 

END !"of procedural” 

PROCEDURAL (D3 = thetar3,Kone,psiqpripri,psidpripri) 
CALL 

simpsond(thetar3,Kone,psiqpripri,psidpripri, D3) 

END !"of procedural" 

PROCEDURAL (Q4 = thetar4,Kone,psiqpripri,psidpripri) 
CALL 

simpsonq (thetar4, Kone, psiqpripri, psidpripri, Q4) 

END !"of procedural" 

PROCEDURAL (D4 = thetar4,Kone,psiqpripri,psidpripri) 
CALL 

simpsond (thetar4, Kone, psiqpripri, psidpripri, D4 ) 

END !"of procedural" 

PROCEDURAL (Q5 = thetarS,Kone,psiqpripri,psidpripri) 
CALL 

simpsonq(thetar5,Kone,psiqpripri,psidpripri, Q5) 

END i"of procedural" 

PROCEDURAL (D5 = thetarS,Kone,psiqpripri,psidpripri) 


CALL 

simpsond(thetar5, Kone, psiqpripri,psidpripri, D5) 

END !"of procedural" 

!"now add all the components together” 
iqsrcom = (u/(PI*4.))*(Ql+4.*Q2+2.*Q3+4.*Q4+Q5) 
idsrcom - (u/(PI*4.))*(Dl+4.*D2+2.*D3+4.*D4+D5) 


average" 


!"add the commutating and conduction components to get 
!"iqsr and idsr" 
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iqsr = iqsrcom+iqsrcond 


idsr = idsrcom+idsrcond 


END !"of derivative" 


END !"of dynamic" 

END !"of program" 

SUBROUTINE simpsonq(thetarsimp,Kone,psiqpripri,psidpripri,iqsrout) 


* subroutine to calculate components for simpson method of determining 

* the q commutating currents 

* 


COMMON / forsub / rad30,rad60,radl20,Xqpripri,Xdpripri,del,PI 

* 

* declare variables 

* 

real ias, parti, part2, part3, Kone, iqsrout, del 

* 

* break the equation for ias down into manageable parts 

* 

partl=Kone-sqrt(3.)*(psiqpripri*cos(thetarsimp+rad30)+psidpripri 
+ *sin(thetarsimp+rad30)) 

* 

part2=(Xqpripri-Xdpripri)*cos(2.*thetarsimp+rad60)+Xqpripri 
+ +Xdpripri 

* 

part3=((Xqpripri-Xdpripri)*cos(2.*thetarsimp+radl20)+0.5 
+ *(Xqpripri+Xdpripri))*dcl 

★ 

* combine the 3 parts to determine ias 

* 

ias=(partl/part2)-(part3/part2) 

* 

* now solve for iqsr at the inputed thetar 

* 

iqsrout=((2.*sqrt(3.))/3.)*(ias*cos(thetarsimp+rad30)-del 
+ *sin(thetarsimp)) 

* 

end 

SUBROUTINE simpsond(thetarsimp,Kone,psiqpripri,psidpripri,idsrout) 
************************************************************************ 
★ 

* subroutine to calculate components for simpson method of determining 

* the d commutating currents 


COMMON / forsub / rad30,rad6Q,radl20,Xqpripri,Xdpripri,del,PI 

* declare variables 

* 

real ias, parti, part2, part3, Kone, idsrout, del 

* 

* break the equation for ias down into manageable parts 
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partl=Kone-sqrt(3.)*(psiqpripri*cos(thetarsimp+rad30)+psidpripri 
+ *sin(thetarsimp+rad30)) 

* 

part2= (Xqpripri-Xdpripri)*cos(2.*thetarsimp+rad60)+Xqpripri 
+ +Xdpripri 

* 

part3=((Xqpripri-Xdpripri)*cos(2.*thetarsimp+radl20)+0.5 
+ *(Xqpripri+Xdpripri))*dcl 

★ 

* combine the 3 parts to determine ias 

* 

ias=(partl/part2)-(part3/part2) 

* 

* now solve for idsr at the inputed thetar 

★ 

idsrout=((2.*sqrt(3.))/3.)*(ias*sin(thetarsimp+rad30)+dcl 
+ *cos(thetarsimp)) 

* 

end 

REAL FUNCTION findu(Kone,psiqpripri,psidpripri) 


* function subprogram to calculate u using the 1 half-it' method 

★ 

COMMON / forsub / rad30,rad60,radl20,Xqpripri,Xdpripri, del, PI 

* 

* declare variables to be used 

★ 

logical condition 

real utemp, fcom, fpricom, Kone, del, ul, u2, u 

* 

ul=0.0 
u2=rad60 
u=0.0 

fcom=Kone-sqrt(3.)*(-psiqpripri*sin(beta+u)+psidpripri 
+ *cos(beta+u))-((Xqpripri-Xdpripri)*cos(2.*beta+2.*u 

+ -radl20)+0.5*(Xqpripri+Xdpripri))*dcl 

10 if (ABS(ul-u2) .gt. 0.001) then 
u = 0.5*(ul+u2) 

fcom=Kone-sqrt(3.)*(-psiqpripri*sin(beta+u)+psidpripri 
+ *cos(beta+u))-((Xqpripri-Xdpripri)*cos(2.*beta+2.*u 

+ -radl20)+0.5*(Xqpripri+Xdpripri))*dcl 

if (fcom .le. 0.00) then 
ul=u 
else 
u2=u 
endif 
goto 10 
endif 
findu=u 
END 

REAL FUNCTION findbeta(psidpripri,psiqpripri, alfahat) 


* function subprogram to calculate beta 
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COMMON / forsub / rad30,rad60,radl20,Xqpripri,Xdpripri, del,PI 


* declare variables to be used 

* 

real betatop, betabot, betamin, alfahat, fbmin 

★ 

betatop = PI 
betabot = -PI 
betamin = 0.0 

* 

* calculate the function to determine if beta is greater than betamin 

* 

fbmin = SQRT(3.)*(psiqpripri*cos(betamin)+psidpripri*sin(betamin)) 
+ +2.*dcl*(Xqpripri-Xdpripri)*sin(2.*betamin-rad60) 

10 if (ABS(betatop-betabot) ,gt. 0.001) then 
betamin = (betatop+betabot)/2. 

fbmin = SQRT(3.)*(psiqpripri*cos(betamin)tpsidpripri 
+ *sin (betamin))+2.*dcl*(Xqpripri-Xdpripri)*sin(2.*betamin 

+ -rad60) 

if (fbmin .le. 0.0) then 
betabot = betamin 
else 

betatop = betamin 
endif 
goto 10 
endif 

★ 

* now determine beta based on the above determination 

* 

findbeta = alfahat + betamin 
END 


B. .CMD FILE 


I"tues.cmd" 

!"Greg Greseth - 15 Feb 1999" 

!"This file plot the output of a synchronous generator, phase- 
controlled" 

!"rectifier model." 


s strplt = .t. 
s calplt = .f. 

s devplt = 1 

s ppoplt = .f. 
s xinspl = 6 
s weditg = .f. 


! "one variable per x-axis" 

! "1 for X-windows" 

! "5 for postscript" 

! "true rotates plot 90 deg" 

! "x-axis plot units" 

! "false suppresses data write" 



s wesitg = .f. 
s nrwitg = .f. 

s alcplt = .f. 
s tstop = 40. 
s maxt=1.0e-4 
s cint=1.0e-3 


! "each time SCHEDULE occurs" 

! "true enables accumulation of data" 
! "after a CONTIN" 

! "false makes all plots dark" 


prepare 

t,Vdc,Te,psikqpri,psikdpri,psifdpri, iqsr,idsr, u, thetar,beta, alpha,phi 

PROCED startup 
s ialg=5 

s devplt = 1 
start 

s pit = 1 

plot Vdc,Te,psikqpri,psikdpri 
s pit = 2 

plot psifdpri,iqsr, idsr,u 
s pit = 3 

plot beta, phi, alpha 

END 

PROCED testalfa 

s tstop = 60 
s ialg=5 
start 

d Vdc,Te,psikqpri,psikdpri, psifdpri, iqsr, idsr, u,beta,phi,alpha 

END 
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APPENDIX B: DERIVATION OF SYNCHRONOUS MACHINE VOLTAGE 

REGULATION 


The reduced-order (stator transients ignored) electrical equations of a 3-phase 
synchrounous machine with a single q and d-axis damper winding may be expressed by: 



(B-l) 


(B-2) 


(B-3) 

^qdr = H * qds + Al qdr 

(B-4) 


where: 




i 



l fd 



(B-5) 


(B-6) 


With the stator currents referenced into the machine and assuming a purely resistive load, 
the terminal voltages may be expressed solely in terms of the terminal currents: 


where: 


~ Rfoad 0 
0 — R-had . 


(B-7) 


(B-8) 


Upon substituting (B-7) into (B-l) and solving for W qds , we get: 

Xds=-B-\A-K)r qds (B-9) 

Further, (B-9) can be substitued into (B-3) so that the stator currents may be solved for in 
terms of the rotor currents: 
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(B-10) 


i qds ^ i qdr 


where: 


l = [-b~\a-k)-f] 1 g 

Equation (B-10) is then plugged into (B-4) to yield: 

T^=(Fr+j)r^ 

and solving for i qdr : 

T qd = (HL+Jy'% dr 

Finally, substituting (B-13) into (B-2), we get: 

^ = c(j7z + J)-% + D^ 

and since = V M = 0, the referred rotor voltage vector may be re-expressed as: 


(B-ll) 

(B-12) 

(B-13) 

(B-14) 


V qdr = NV fd = 


fd 


(B-15) 


The output of the machine is the terminal voltages. Equations (B-7), (B-10), and 
(B-13) may be combined to express the terminal voltages in terms of the rotor states: 

V^=KL(HL+Jy l % dr (B- 16 ) 

Once we assume that the turbine governor is holding co r constant, equations (B-14) and 
(B-16) become a set of linear ordinary differential equations. However, the peak value of 
the terminal voltages is related to V qs and by: 

^,M=V(^) 2 +(^) 2 (B-17) 


Therefore, to consider a single-input ( V ' Jd ), single-output (V qds pk ) system, the output 
relationship must be linearized about an operating point: 

=MAV q y (B-18) 

where: 
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f dso _ 

r \2 , /T/r \2 


yKoY+iVLY yl^soY+iVL) 

Finally, the state-space model of the generator is given by: 

= C 2 A^ 

where: 


(B-19) 


■A% dr =A 2 A% dr+ B 2 AV fd 

(B-20) 

= C 2 A% r 

(B-21) 

a 2 =-n-'c(HL+Jy 1 

(B-22) 

B 2 = D~ X N 

(B-23) 

C 2 = MKL(H L + J)~ x 

(B-24) 


A final point lies in determining the operating point values for (B-19). This can 
be accomplished assuming that for rated load, V /d has been selected to provide the rated 
terminal voltage. In the steady state: 


v T /e --A R V 

*qdro D 1 V Jdo 

Substituting (B-25) into the output equation (B-16), we find: 

v; dS 0 =-KL(HL+Jy'A 2 -%v /c 

or: 


where A 3 is a 2-by-l vector: 


V = A V' 

V qdso Sd 


7 _ 40 ) 

Ay — — 

4C 2) 


(B-25) 


(B-26) 


(B-27) 


(B-28) 


Thus, from (B-17) and (B-27): 


V q *, pk =ylM'Y+M2yV /d (B-29) 

Therefore, for a given load resistance ( K ), the appropriate field voltage to give the rated 
terminal voltage is found from: 
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1 


(B-30) 




2? 


F r 

qds,pk,rat 


This value may then be substitued back into (B-27) to yield the desired operating point 
values. 


At this point, the design of the field regulator can commence. To illustrate the 
need for compensation, the open-loop eigenvalues of A 2 from Equation B-20 can be 
calculated assuming a load resistance of 1.0 pu. These are found to be: 


Xy =-19.71 
Z 2 = -198.61 


X 3 = -0.54 


The A 3 eigenvalue introduces a fairly long recovery time. In addition, open-loop, 
the operator would need to know the load in order to apply the correct V fd to maintain the 

desired pk . For instance, with R load = 1.0, the required field voltage to yield pk = 
1.0 is: 


V' fd = 0.0012789 
or 


4=2.039=^?% (B-31) 

r* 

As the load power is decreased, this excitation voltage must be lowered to maintain 
Vqds, P k = 1-0. For the above condition, the required stator voltage operating points are 
found to be: 

V; o = 0.6813 pu 
VL = 0.7320 pu 

The system can be linearized about these operating points and a transfer function derived 
between &V' fd and AV^ pk : 
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(B-32) 


G(s) = 


= 219.7(5 + 271)(s + 27.7) 
AV' fd (s +19.7)0 + 199)(s + 0.54) 


The desired control loop is then illustrated in Figure B-l. The circuit can be 
analyzed to yield the closed-loop transfer function: 


and the steady-state error: 


= n c n g 

W q *, P k n c n g +d c d g 


E = w£l 

n c n g +d c d g s=Q 


where n c , n g , d c and d g are polynomials in s. 


(B-33) 


(B-34) 


+ 


1 

< 


d c c 2 s +1 


1 

D 


Figure B-l Field Regulator Control Loop 


The control gains c, and c 2 must be selected to yield satisfactory steady-state 
error and acceptable closed-loop pole locations. Plugging in the transfer function 
elements into Equation B-34 yields: 

1 


£ = 


l + 799c, 


(B-35) 


If the output voltage is to be within 0.5% of the commanded value, this imposes the 
constraint: 

c,> 0.223 (B-36) 
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If it is assumed that c x = 0.25 for convenience, c 2 may be selected based on the closed- 
loop poles of Equation B-33. Note, the zeros of Equation B-33 are fixed by the zeros of 
n g and once c, is fixed, the steady-state error is fixed. Table B-l illustrates how the 
closed-loop poles vary as c 2 is changed. 


C 2 

Pole 1 

Pole 2 

Poles 3 & 4 

1000 

-199 

-19.7 

-0.2697j0.18 

10 

-199 

-19.9 

-0.2357j3.23 

1 

-199 

-21.0 

-0.0847j 9.97 

0.1 

-200 

-25.1 

-2.067j28.7 

0.01 

-212 

-30.6 

-38.237j70.23 

0.002 

-457 

-32.6 

-1147j28.3 

0.001 

-948 

-33 

-89.1 & -149 


Table B-l Closed-loop Poles Versus Gain c 2 


As the gain is decreased, the poles 3 & 4 move from slow, lightly damped 
locations to positions comfortably in the left-half plane. An attractive compensator 
transfer function is given by: 


or 


AV, 


fd 




0.25 

0.0025 + 1 


(B-37) 


Ae 


xfd 


400 


wi Pk 


0.0025 + 1 

One final point to understand is that this design is valid for one operating 


(B-38) 


condition ( R load — 1.0 pu). As illustrated by Table B-2, the steady-state value of e x j d 
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must be varied to maintain rated output voltage. The corresponding transfer functions 
can be evaluated and analyzed for suitability. 


f 

I 

e x jd (P u ) 

1.0 

2.039 

1.2 

1.771 

1.4 

1.601 

1.6 

1.478 

1.8 

1.388 

2.0 

1.321 

5.0 

1.057 

10.0 

1.015 


Table B-2 Excitation Voltage Required to Yield Rated Output Voltage 
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APPENDIX C: SIMPLIFIED DC LINK ACSL MODEL 


A. NON-LINEAR .CSL FILE 


PROGRAM !"simplenonlin.csl" 

!"Greg Greseth — 8 April 1999" 

!"This file contains the non-linear system equations for a simplified" 
!"DC link model. The simulated event is a 10%-100%-10% power 
transient." 


INITIAL 

ALGORITHM ialg = 5 
MAXTERVAL maxt = l.e-5 
MINTERVAL mint = l.e-6 
CINTERVAL cint = 5.e-3 
CONSTANT Le = 0.0005 
CONSTANT Ce = 0.008265 
CONSTANT Re = 0.01 
CONSTANT tstop =5.0 
CONSTANT P = 20000. 

CONSTANT Ves = 1102. 

END !"of initial 

DYNAMIC 

TERMT(t.ge.tstop, 'Termination on time limit') 

DERIVATIVE 

P_star = P + 180000.*STEP(2.0) - 180000.*STEP(3.5) 
ir = INTEG(-(Vdci/Le)-(Re*ir)/Le+Ves/Le,181.82) 
Vdci = INTEG(ir/Ce-(P_star)/(Ce*Vdci),1100.) 

END !"of derivative" 

END !"of dynamic" 

END !"of program" 
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B. 


NON-LINEAR .CMD FILE 


! "simplenonlin.cmd" 

!"Greg Greseth - 8 April 1999" 


! "This file plots the inductor current and capacitor voltage" 
! "for the non-linear simplified DC link model" 


s strplt = .t. 
s calplt = . f. 

s devplt = 1 

s ppoplt = .f. 
s xinspl = 6 
s weditg = .f. 
s wesitg - . f. 
s nrwitg = .f. 

s alcplt = .f. 


! "one variable per x-axis" 


! "1 for X-windows" 

! "5 for postscript" 

! "true rotates plot 90 deg" 

! "x-axis plot units" 

! "false suppresses data write" 
! "each time SCHEDULE occurs" 


! "true enables accumulation of data" 
! "after a CONTIN" 

! "false makes all plots dark" 


prepare t,ir,Vdci 


PROCED nonlinear 
start 

plot ir,Vdci 

END 


C. LINEAR .CSL FILE 


PROGRAM !"simplelin.csl" 

!"Greg Greseth - 8 April 1999" 

!"This file contains the linearized system equations for a simplified" 
!"DC link model. The simulated event is a 10%-100%-10% power 
transient." 

INITIAL 

ALGORITHM ialg = 5 
MAXTERVAL maxt = l.e-5 
MINTERVAL mint = l.e-6 
CINTERVAL cint - 5.e-3 
CONSTANT Le = 0.0005 
CONSTANT Ce = 0.008265 
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CONSTANT Re = 0.01 
CONSTANT tstop =1.5 
CONSTANT P = 20000. 

CONSTANT Vdci = 1100. 

CONSTANT delta_Ves =5. 

END !"of initial 

DYNAMIC 

TERMT(t.ge.tstop, 'Termination on time limit') 

DERIVATIVE 

P_star = P + 180000.*STEP(0.5) - 180000.*STEP(3.0) 

delta_ir = INTEG<-(delta_Vdci/Le)-(Re*delta_ir)/Le+delta_Ves/Le,0) 

delta_Vdci = INTEG(delta_ir/Ce+(P_star*delta_Vdci)/(Ce*Vdci**2),5) 

END. !"of derivative" 

END !"of dynamic" 

END !"of program" 

D. LINEAR .CMD FILE 


!"simplelin.cmd" 

I"Greg Greseth - 8 April 1999" 

! "This file plots the linearized inductor current and capacitor 1 ' 
•"voltage for the simplified DC link model" 


s strplt = .t. 
s calplt - .f. 

s devplt = 1 


s ppoplt = . f. 
s xinspl = 6 
s weditg = .f. 
s wesitg = . f. 
s nrwitg = .f. 


! "one variable per x-axis" 

! "1 for X-windows" 

! "5 for postscript" 

! "true rotates plot 90 deg" 

! "x-axis plot units" 

! "false suppresses data write" 
! "each time SCHEDULE occurs" 


! "true enables accumulation of data" 
! "after a CONTIN" 


s alcplt = .f. 


"false makes all plots dark" 
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prepare t, delta_ir, delta_Vdci 


PROCED linear 
start 

plot delta_ir,deltaJVdci 

END 
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APPENDIX D: FILTER-CONVERTER ACSL MODEL 


A. .CSL FILE 


!"hardsw.csl" 

!"written by John Ciezki" 

!"modified by Greg Greseth - 29 April 1999" 

’"last modified - 1 May 1999" 

!"This models a DC-DC converter (buck chopper) with an DC link filter 
at" 

!"the input." 


PROGRAM 


INITIAL 


MAXTERVAL maxt = 2.0e-7 
MINTERVAL mint * 1.0e-7 


! "maximum integration step size" 


CINTERVAL cint = 
ALGORITHM ialg = 

NSTEPS nstp = 
CONSTANT tstop 


l.Oe-3 !"data communication interval" 

5 !"integration algorithm" 

! "4—R.K. 2nd, 5—R.K. '4th" 

1 

= 0.25 !"stop point for Integration" 


!"Buck Parameters" 
CONSTANT Lbuck - 1.0e-3 
CONSTANT Riser =0.0 
CONSTANT Vdidrop =0.0 
CONSTANT Vswdrop = 0.0 
CONSTANT Cbuck = 2.0e-3 
CONSTANT Tramp = 250.0e-6 


"main inductance" 

"Lbuck series resistance" 
"diode forw conduction drop" 
"switch forward drop" 

"main output capacitance" 
"switching period (4kHz)" 


!"-define and initialize the main switch conduction status" 

LOGICAL sw 

sw = .false. !"T sw conducts, F sw blocks" 

i»-define and initialize variable which keeps track of" 

!" whether the converter is in contin or discontin mode" 

LOGICAL contmode 
contmode = .true. 


’"link filter parameters" 
CONSTANT Cf = 8265.Oe-6 
CONSTANT Lf = 0.5e-3 
CONSTANT Rs = 0.01 


’"filter capacitance" 
! "filter inductance" 

! "source resistance" 


!"state variable initial conditions" 

CONSTANT iLic = 0.0 !"Buck inductor current" 
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CONSTANT Voic =0.0 !"Buck output voltage" 

CONSTANT Vcfic = 1100.0 !"filter output voltage" 

CONSTANT iLfic = 0.0 !"filter inductor current" 

!"Load Parameters" 

CONSTANT Rmax = 4.5125 ! "100% load" 

CONSTANT Rmin = 45.125 !"min load" 

!"start-up" 

CONSTANT trampup = 0.1 ! "time to ramp from 0 to vrefmax" 

CONSTANT vrefmax = 950.0 !"max ref output voltage" 

LOGICAL rampon 
rampon = .TRUE. 

!"define the input voltage" 

CONSTANT Vin =1100.0 

END !"of initial" 

DYNAMIC 

TERMT (t .GE. (tstop-0.5*cint)) 

DERIVATIVE 

!"Set the load resistance" 

Rload = Rmax 

SCHEDULE rampoff .XP. t-trampup 
IF (rampon) THEN 

Vref = Vrefmax*t/trampup 

ELSE 

Vref = Vrefmax 

ENDIF 

‘"BUCK MODELING EQUATIONS" 

11 -determine conduction status of swl" 

PROCEDURAL(sw,isw=reframp,duty,iL) 

IF (duty .LE. reframp) THEN 

sw = .false. !"switch open" 

isw = 0.0 !"source current zero" 

ELSE 

sw = .true. !"switch closed" 
isw = iL 

ENDIF 

END !"of procedural" 

" given conduction status of sw and whether current is" 

" continuous or not, establish derivative of iL" 

PROCEDURAL(piL=sw,contmode,Vcf,Vo,iL,isw) 

• ” if sw gated, assume can conduct pos current" 

IF (sw) THEN 

contmode = .true. 

piL = (-Rlser*iL+Vcf-Vswdrop-Vo)/Lbuck 
ELSE 

"-if iL GT zero, diode conducts" 

IF (contmode) THEN 
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piL = (-Rlser*iL-Vdidrop-Vo)/Lbuck 
ELSE 

piL = 0.0 !"discontinuous mode" 

ENDIF 

ENDIF 

END !"of procedural" 

!"link filter state variables integrated" 

Vcf = INTEG((iLf-isw)/Cf, Vcfic) 

iLf = INTEG((Vin-Vcf-(Rs*iLf))/Lf, iLfic) 

!"establish the derivative of Buck output voltage" 
pVo = (iL - io)/Cbuck 

!"integrate the state variables" 

iLub = INTEG(piL, iLic) 

iL = BOUND(-1.Oe-6, 1.0e6, iLub) 

io = Vo/Rload 

Vo = INTEG(pVo, Voic) 

"-form the reference ramp waveform (0 to 1)" 

tt = mod(t. Tramp) 
reframp = tt/Tramp 

"-schedule discontinuous mode when iL tries go neg" 

SCHEDULE DCM ,XN. iL 

!"Control Algorithm for Buck" 

CONSTANT hv = 0.00914 
!CONSTANT hi = -0.003299 
CONSTANT hi = -0.003208 
CONSTANT hn = 6.093 

CONSTANT nic = 0.0 ' 

CONSTANT intbnd = 1000.0 
CONSTANT minduty =0.0 
CONSTANT maxduty =0.99 
CONSTANT Kscv = 1.000 
CONSTANT Ksci = 1.000 
va = Vref-Vo 
ioave = Ksci*io 
iLave = Ksci*iL 
ia = iLave-ioave 

pn = va 

nub = INTEG(pn, nic) 
n = BOUND(-intbnd, intbnd, nub) 
dss = Vref/Vin 

dub =dss + Kscv*hv*va + Ksci*hi*ia + Kscv*hn*n 
duty = BOUND(minduty, maxduty, dub) 

END ! "of derivative" 

vout = Rload*io 

"-establish transition to discontinuous mode Buck" 

DISCRETE DCM 

contmode = .false. 


{"integral of volt err" 

{"scale factor for voltage sensor" 
{"scale factor for current sensor" 
I "proportional voltage term" 

I "proportional current term" 

{ "integral voltage term" 
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END 


! "of dcm 


■ DISCRETE rampoff 

rampon = .FALSE. 

END 

END ! "of dynamic" 

END ! "of program" 

B. .CMD FILE 


!"hardsw.cmd" 

!"Greg Greseth - 1 May 1999" 

!"This file plots the output for a filter-converter model" 


s strplt = .t. 
s calplt = .f. 


! "one variable per x-axis" 


s devplt = 1 


s ppoplt = .f. 
s xinspl = 6 
s weditg = ,f. 
s wesitg = .f. 


! "1 for X-windows" 

! "5 for postscript" 

I "true rotates plot 90 deg" 

! "x-axis plot units" 

! "false suppresses data write" 
! "each time SCHEDULE occurs" 


s nrwitg = .f. 


! "true enables accumulation of data" 
i "after a CONTIN" 


s alcplt = .f. 


'false makes all plots dark" 


prepare t,Vcf,iLf,iL,io,Vo,duty 

PROCED startup 
s ialg=5 
s maxt=1.0e-6 
s cint=5.Oe-4 
start 

s pit = 1 
plot Vcf,iLf,duty 
s pit = 2 
plot iL,io,Vo 

END 
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APPENDIX E: MATLAB CODE FOR IMPEDANCE CALCULATIONS 


A. HARD SWITCHING MODEL 


% findxferhard.m 
% 

% Greg .Greseth - 30 March 99 
% modified 19 April 99 

% 

% This program will find the transfer function for a linearized buck 
model 

% in order to determine T2 and T3 for the hard switching case. For 
analysis 

% purposes, delta Vref is assumed to be 0. 

% 


% define the constants 
Pin - 200e3; 

Vrefo =950; % =Vout 

Vino = 1100; 
do = 0.8636; 

L = le-3; 

C = 2e-3; 

R = 45.125; % =Rcrit 

initdeltaVin = 5; 

K1 = -Vrefo/(Vino) A 2; 
ilo = Pin/Vino; 

% characteristic equation coeffiecients 

d2 = 3639.5; 

dl = 5526354.8; 

d_zero = 3351129059.4; 

% calculate the gains 

hv = {(L*C)/Vino)*(dl-(l/(L*C))); %hv = 0.00914; 

hn = (d__zero*L*C) /Vino; %hn = 6.093; 

hi = -(L/Vino)*(d2-(1/(R*C))); %hi = 0.003299; 

Tlnum = [hi*C, hv, hn]; 

% Define the A matrix 

A = [(Vino*hi)/L,(-1/L)-(Vino^hv)/L-(Vino^hi)/(R*L),(Vino*hn)/L;... 
1/C,-1/(R*C),0;... 

0 ,- 1 , 0 ]; 

% Define the B matrix 
B = [-Vrefo/(L*Vino)+do/L;0;0]; 
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% Define the D matrix, which is always 0 
D = 0; 

% Define C matrix in order to find T2 
C_T2 = [1,0,0]; 

% Find the transfer function which will define T2 
[T2num,T2den] = ss2tf(A,B,C_T2,D) 

I Redefine C matrix in order to find T3 
C_T3 = [0,1,0]; 

% Find the transfer function which will define T3 
[T3num,T3den] = ss2tf(A,B,C_T3,D) 

% Determine the Zin numerator 
Zin_num_hard = [T2den, 0] 

% Determine the parts of the Zin denomenator 

parti = [do*T2num, 0]; 

part2 = [ilo*Kl*T3den, 0]; 

part3 = ilo*conv(Tlnum, T3num); 

part3 = part3(1,2:6); 

Zin_den_hard = partl+part2+part3 

B. SOFT SWITCHING MODEL 

% findxfersoft.m 
% 

% Greg Greseth - 30 March 99 
% modified 19 April 99 

% 

% This program will find the transfer function for a linearized buck 
model 

% in order to determine T2 and T3 for the soft switching case. For 
analysis 

% purposes, delta Vref is assumed to be 0. 

% 

% define the constants 
Pin = 200e3; 
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% =Vout 


Vrefo = 950; 

Vino = 1100; 
do = 0.8636; 

L = 200e-6; 

C = 400e-6; 

R = 45.125; % =Rcrit 

initdeltaVin = 5; 

K1 = -Vrefo/(Vino)^2; 
ilo = Pin/Vino; 


% characteristic equation coefficients 

d2 = 7299; 

dl = 22194879.2; 

d_zero = 27000116071.2; 

% calculate the gains 

hv = ((L*C)/Vino)*(dl-(l/(L*C))); %hv = 0.00070508; 

hn = (d_zero*L*C)/Vino; %hn = 1.9636; 

hi = -(L/Vino)*(d2-(1/(R*C))); %hi = 0.001317; 

Tlnum = [hi*C, hv, hn]; 

% Define the A matrix 

A = [(Vino*hi)/L,(-1/L)-(Vino*hv)/L-(Vino*hi)/(R*L),(Vino*hn)/L 
1/C,-1/(R*C),0;... 

0 ,- 1 , 0 ]; 

% Define the B matrix 
B = [-Vrefo/(L*Vino)+do/L;0;0]; 

% Define the D matrix, which is always 0 
D = 0; 


% Define C matrix in order to find T2 
C_T2 = [1,0,0]; 

% Find the transfer function which will define T2 
[T2num,T2den] = ss2tf(A, B, C_T2,D) 

% Redefine C matrix in order to find T3 
C_T3 = [0,1,0]; 

% Find the transfer function which will define T3 
[T3num,T3den] = ss2tf(A,B,C_T3,D) 

% Determine the Zin numerator 
Zin num soft = [T2den, 0] 
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% Determine the parts of the Zin denomenator 

parti = [do*T2num, 0]; 

part2 = [ilo*Kl*T3den, 0] ; 

part3 = ilo*conv(Tlnum, T3num); 

part3 = part3(1,2:6); 

Zin_den_soft = partl+part2+part3 


C. NYQUIST PLOT FILE 


% stability.m 
% 

% Greg Greseth - 19 April 99 

% 

% This program will find the Nyquist plot of the impedance ratio 
% between the filter output impedance and the converter input 
% impedance. The desired result is to have the magnitude of the 
% ratio be less than one over all range of angular frequency. 

% Define the functions being analyzed: 

% Zs is output impedance of the filter 
% Zi is input impedance of the converter 

% hard indicates values used for hard switching converter 
*6 soft indicates values used for soft switching converter 

% IMPORTANT !!!!!!!!!!! 

% Both findxferhard.m and findxfersoft.m must be run with the desired 
% converter values must be run first to ensure that the input 
% impedance vectors are in the matlab variable list 

% define the filter parameters and source resistance 

Rs = 0.01; % Source resistance 

Cf - 6000e-6; % Filter capacitance - 8265e-6 

Lf = 0.5e-3; % Filter inductance - 0.5e-3 

% calculate the filter output impedance 

Zsnum = [1/Cf, Rs/(Cf*Lf)]; % numerator 

Zsden = [1, Rs/Lf, 1/(Cf*Lf)]; % denomenator 

% calculate the hard switching ratio for Nyquist test 

hardnum = conv(Zsnum,Zin_den_hard); 
hardden = conv(Zsden,Zin_num_hard) ; 

% calculate the soft switching ratio for Nyquist test 
softnum = conv(Zsnum,Zin_den soft); 
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softden = conv (Zsden, Zin_num_soft) ; 

% Nyquist plot for hard switching case 
figure (1) 

nyquist(hardnum,hardden) 

% Nyquist plot for soft switching case 
figure (2) 

nyquist(softnum,softden) 
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APPENDIX F: ACSL MODEL FOR LYAPUNOV-BASED STABILITY 


A. .CSL FILE 


PROGRAM !"oneload.csl" 

!"Greg Greseth - 13 May 1999" 

!"This file finds the load current and load voltage for a single load 
system” 

!”for the Lyapunov-based stability analysis. The non-linear system 
equations” 

!”for the simplified model are utilized. The simulation event models a” 
!"disturbance (drop) in the source voltage." 

INITIAL 

CONSTANT R=0.01 
CONSTANT L1=0.0005 
CONSTANT Cl=.008265 
CONSTANT Pl=200000.0 
ALGORITHM ialg = 5 
CONSTANT tstop =1.0 
NSTEPS nstp=l 
MAXTERVAL maxt = l.Oe-4 
MINTERVAL mint = 1.0e-7 
CINTERVAL cint = 4.0e-3 
CONSTANT DROP =55.0 
CONSTANT VDCSS = 1100.0 

END !"of initial" 

DYNAMIC 

DERIVATIVE 

VDC=VDCSS-DROP*step(0.5) !"supply voltage" 

!"initial conditions" 

CONSTANT I1IC=190.0 

CONSTANT V1IC=1000.0 


"resistance value" 

"inductor 1 value" 

"capacitor 1 value" 

"load 1 power" 

"Runge-Kutta 4th order" 

"stop time for integration" 

"max integration step size" 
"min integration step size" 
"data communication interval" 
"source voltage disturbance" 
"source voltage steady state" 


! "current" 


I1=INTEG(-R/L1*(II)-V1/L1+VDC/L1,11IC) 

! "conditions for VI" 

IF(VI.LE.0.01) THEN 
pVl=Il/Cl 

ELSE 
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PV1=I1/C1-P1/(C1*V1) 


ENDIF 

!"voltages" 

Vl=INTEG(pVl,VlIC) 

TERMT(T . GE .TSTOP, ’Termination on time limit 1 ) 
END !"of derivative" 

END !"of dynamic" 

END !"of program" 


B. .CMD FILE 


!"oneload.cmd" 

!"Greg Greseth - 13 May 1999" 

!"This file plots the load current and load voltage for a single load 
system" 

!"for the Lyapunov-based stability analysis" 

s strplt = *t. !"one variable per x-axis" 

s calplt — .f. 

s devplt - 1 !"1 for x-windows" 

!"5 for postscript" 

prepare t, II, VI 
procedure plotl 
start 

plot II /xtag=*time* /tag=*Load 1 Current’ VI /tag=*Load 1 
Voltage* 

end 
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APPENDIX G: ACSL MODEL FOR INTEGRATED SYSTEM 


A. .CSL FILE 


PROGRAM !"system.csl" 

!"Greg Greseth, LT, USN - 13 May 99" 

!"Last Revised - 26 May 99 

!"This ACSL file simulates a reduced order model for a 

synchronous" 

!"machine and three phase rectifier, with the output through 

a DC" 

!"link filter and a buck chopper with a resistive load." 

!"Necessary inputs/known values are:" 

!"Output from the model is:" 


INITIAL !"Section in which constants and program 

parameters are defined" 



ALGORITHM 

ialg = 5 

!"Runge-Kutta 4th order 


CONSTANT tstop =0.25 

!”stop time for 

integration 

tr 

MAXTERVAL 

maxt = 5.e-7 

!"max integration step 

size” 


MINTERVAL 

mint = l.e-6 

!”min integration step 

size” 


CINTERVAL 

cint = 5.e-4 

!"data communication 

interval” 


CONSTANT wb = 377. 

CONSTANT rkqpri = 0.0613 
CONSTANT rkdpri =0.24 
CONSTANT rfdpri = 0.00111 
CONSTANT J = 1. 

CONSTANT T1 = 1. 

CONSTANT B = 0. 

CONSTANT P = 2. !"number of poles" 

CONSTANT Xls =0.08 

CONSTANT Xmq = 1.00 

CONSTANT Xlkqpri = 0.330 

CONSTANT Xmd =1.77 

CONSTANT Xlfdpri = 0.137 

CONSTANT Xlkdpri = 0.334 

CONSTANT rs = 0.00515 

CONSTANT exfdo = 0.0042 

CONSTANT alfahat =0.0 
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CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 


PI = 3.141593 
radl20 = 2.094395 
rad60 - 1.047198 
rad30 = 0.523599 
Vpkstar = 1.654 
confact = 400.0 
sterm = 0.002 


call" 


u = 0.05 !"initialize u for use in first findu function 


{"define the subtransient reactances" 

Xqpripri = Xls+ (Xmq*Xlkqpri) / (Xlkqpri+Xrriq) 
Xdpripri = 

Xls+ (Xmd*Xlfdpri*Xlkdpri) / (Xmd*Xlfdpri+Xmd*Xlkdpri& 
+X1kdpri*X1fdpri) 


{"Buck Parameters" 
CONSTANT Lbuck = 1.0e-3 
CONSTANT Riser = 0.0 
CONSTANT Vdidrop =0.0 
CONSTANT Vswdrop =0.0 
CONSTANT Cbuck = 2.Oe-3 
CONSTANT Tramp = 250.0e-6 


"main inductance" 

"Lbuck series resistance" 
"diode forw conduction drop" 
"switch forward drop" 

"main output capacitance" 
"switching period (4kHz)" 


* -define and initialize the main switch conduction status" 

LOGICAL sw 

sw = .false. !"T sw conducts, F sw blocks" 

i«-define and initialize variable which keeps track of" 

I" whether the converter is in contin or discontin mode" 

LOGICAL contmode 
contmode = .true. 

{"link filter parameters" 

CONSTANT Cf = 8265.0e-fe {"filter capacitance" 

CONSTANT Lf = 0.5e-3 {"filter inductance" 

CONSTANT Rsrc = 0.01 {"source resistance" 


{"state variable initial conditions" 

CONSTANT iLic = 0.0 {"Buck inductor current" 

CONSTANT Voic =0.0 {"Buck output voltage" 

CONSTANT Vcfic = 1100.0 {"filter output voltage" 
CONSTANT iLfic = 0.0 {"filter inductor current" 


{"Load Parameters" 

CONSTANT Rmax = 4.5125 {"100% load" 

CONSTANT Rmin = 45.125 {"min load" 

{"start-up" 

CONSTANT trampup 
CONSTANT vrefmax 
LOGICAL ramp on 
rampon = .TRUE. 

{"scale between per unit values and actual values" 


0.1 {"time to ramp from 0 to vrefmax' 

950.0 {"max ref output voltage" 
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CONSTANT currscal = 202.0 
CONSTANT voltseal = 411.0 


COMMON / forsub / 

rad30, rad60,radl20,Xqpripri,Xdpripri, PI,currscal 


END !"of initial section" 

DYNAMIC 

TERMT(t.ge.(tstop-0.5*cint), 'Termination on time limit') 
DERIVATIVE 

!"equations included for field voltage control of the 

generator" 

Vfdpri = (exfd*rfdpri)/Xmd 
Vpk = sqrt(Vqsr*Vqsr+Vdsr*Vdsr) 

Vdiff = Vpkstar - Vpk 

exfdub = confact*realpl{sterm,Vdiff,exfdo) 
exfd = BOUND (-7.0, 7.0, exfdub) 

!"state equations for psikqpri, psikdpri, psifdpri, wr 

and thetar" 

psikqpri - INTEG(-wb*rkqpri*ikqpri, -0.120) 

psikdpri = INTEG(-wb*rkdpri*ikdpri, 1.647) 

psifdpri = INTEG(-wb*rfdpri*ifdpri+wb*Vfdpri, 1.777) 

!"wrm = INTEG(((l./J)*(Te-Tl-B*wrm)), 0.)" 

!"relationship between wrm and state wr 

!"wr = (P/2.)*wrm" 
wr = 377. 

thetar = INTEG(wr, 0.) 

!"thetar = mod(theta,rad60) + rad60d" 

{"electrical torque equation" 

!"Te = ( (3./P)/(4.*wb))*(psidsr*iqsr-psiqsr*idsr)" 

Te = psidsr*iqsr-psiqsr*idsr 

{"solve first for ikdpri, ifdpri and ikqpri" 

ikdpri = (1./((Xlkdpri+Xmd)*(Xlfdpri+Xmd)- 
Xmd**2 .) ) * ( (Xlfdpri+Xmd) & 
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*(psikdpri-Xmd*idsr)-Xmd*(psifdpri-Xmd*idsr)) 

ifdpri = (1./((Xlkdpri+Xmd)*(Xlfdpri+Xmd)-Xmd**2.))*(- 

Xmd*(psikdpriS 

-Xmd*idsr)+(Xlkdpri+Xmd)*(psifdpri-Xmd*idsr)) 

ikqpri = psikqpri/(Xlkqpri+Xmq)- 
(Xmq*iqsr)/(Xlkqpri+Xmq) 

("equations to eliminate psiqsr and psidsr" 

psiqsr = Xqpripri*iqsr+psiqpripri 

psidsr = Xdpripri*idsr+psidpripri 

("define the subtransient flux linkages" 

psiqpripri = (Xmq/(Xmq+Xlkqpri))*psikqpri 

psidpripri = 

((Xmd*psikdpri)/Xlkdpri+(Xmd*psifdpri)/Xlfdpri)& 
/(l.+Xmd/Xlfdpri+Xmd/Xlkdpri) 

("define Eqpripri and Edpripri" 

Eqpripri = (wr*psidpripri)/wb 

Edpripri = -(wr*psiqpripri)/wb 

!"define phi to use to determine alpha" 

phi = ATAN2(Edpripri,Eqpripri) 

("call function to determine beta" 

beta = 

findbeta(psidpripri,psiqpripri,iLf, alfahat) 

!"define alpha" 
alpha = beta - phi 

!"desired reduced order voltage-behind-reactance 

equations" 

Vqsr = rs*iqsr+(wr/wb)*Xdpripri*idsr+Eqpripri 

Vdsr = rs*idsr-(wr/wb)*Xqpripri*iqsr+Edpripri 

!"equation for average DC voltage out of rectifier" 

Vdc = (3.*sqrt(3.)/PI)*(wr/wb)*(sqrt(psiqpripri**2& 
+psidpripri**2))*(cos(alpha))-(3./PI)*(wr/wb)*({(Xqpripri+Xdpripri)/2.)& 
*(iLf/currscal)+(Xdpripri-Xqpripri)*(iLf/currscal)*sin(2.*beta+rad30)) 

("need to define K1 for use in subroutines" 
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Kone = sqrt(3.)*{- 

psiqpripri*sin(beta)+psidpripri*cos(beta))& 

+ { (Xdpripri-Xqpripri)*cos(2.*beta+radl20)-0.5*(Xqpripri+Xdpripri))& 
*(iLf/currscal) 


!"call function to determine u using 'half-it* method" 
u « findu(Kone,psiqpripri,psidpripri) 

!"now with the new u, we can determine conduction 

currents" 


iqsrcond = ((2.*sqrt(3.))/PI)*(iLf/currscal)& 
*(cos(beta+radl20)-cos(beta+u+rad60)) 

idsrcond = ((2.*sqrt(3.))/PI)*(iLf/currscal)& 
*(sin(beta+radl20)-sin(beta+u+rad60)) 


method" 


! "calculate the commutating currents using Simpson 


used" 


!"first, define the 5 different values of thetar to be 


thetarl = beta+rad60 
thetar2 = beta+rad60+u/4. 
thetar3 = beta+rad60+u/2. 
thetar4 = beta+rad60+(3.*u)/4. 
thetarS = beta+rad60+u 

•"call each q and d subroutine 5 times to calculate 

components" 


PROCEDURAL (Q1 = 

thetarl,Kone,psiqpripri,psidpripri, iLf) 


CALL 

simpsonq(thetarl,Kone,psiqpripri,psidpripri,iLf, Q1) 

END !"of procedural" 

PROCEDURAL (D1 = 

thetarl,Kone,psiqpripri,psidpripri, iLf) 

CALL 

simpsond(thetarl,Kone,psiqpripri,psidpripri,iLf, Dl) 

END !"of procedural" 

PROCEDURAL (Q2 = 

thetar2,Kone,psiqpripri,psidpripri, iLf) 
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CALL 

siinpsonq {thetar2, Kone, psiqpripri, psidpripri, iLf, Q2) 
END ! "of procedural 11 
PROCEDURAL (D2 = 

thetar2,Kone,psiqpripri,psidpripri,iLf) 

CALL 

simpsond(thetar2,Kone,psiqpripri,psidpripri,iLf,D2) 
END !"of procedural” 

PROCEDURAL (Q3 = 

thetar3,Kone,psiqpripri,psidpripri,iLf) 

CALL 

simpsonq(thetar3,Kone,psiqpripri,psidpripri,iLf,Q3) 
END !"of procedural" 

PROCEDURAL (D3 = 

thetar3,Kone,psiqpripri,psidpripri,iLf) 

CALL 

simpsond(thetar3,Kone,psiqpripri,psidpripri,iLf,D3) 
END !"of procedural" 

PROCEDURAL (Q4 = 

thetar4,Kone,psiqpripri,psidpripri,iLf) 

CALL 

siinpsonq (thetar4, Kone, psiqpripri,psidpripri, iLf, Q4) 
END !"of procedural" 

PROCEDURAL (D4 « 

thetar4,Kone,psiqpripri,psidpripri,iLf) 

CALL 

simpsond(thetar4,Kone,psiqpripri,psidpripri,iLf,D4) 
END !"of procedural" 

PROCEDURAL (Q5 = 

thetarS,Kone,psiqpripri,psidpripri,iLf) 

CALL 

simpsonq(thetar5,Kone,psiqpripri,psidpripri,iLf,Q5) 
END !"of procedural" 

PROCEDURAL (D5 = 

thetar5,Kone,psiqpripri,psidpripri,iLf) 
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CALL 

simpsond(thetar5,Kone,psiqpripri,psidpripri,iLf,D5) 

END !"of procedural" 

!"now add all the components together" 
iqsrcom = (u/(PI*4.))*(Ql+4.*Q2+2.*Q3+4.*Q4+Q5) 
idsrcom = (u/ (PI*4.))*(Dl+4.*D2+2.*D3+4.*D4+D5) 

! "add the commutating and conduction components to get 

average" 

!"iqsr and idsr" 

iqsr = iqsrcom+iqsrcond 

idsr = idsrcom+idsrcond 

!"Set the load resistance" 

Rload = Rmax 


SCHEDULE rampoff .XP. t-trampup 
IF (rampon) THEN 

Vref = Vrefmax*t/trampup 

ELSE 

Vref = Vrefmax 

END IF 

!"BUCK MODELING EQUATIONS" 

"-determine conduction status of swl" 

PROCEDURAL(sw,isw=reframp,duty,iL) 

IF (duty .LE. reframp) THEN 

sw = .false. !"switch open" 

isw = 0.0 '"source current zero" 

ELSE 

sw = .true. '"switch closed" 
isw = iL 

ENDIF 

END !"of procedural" 

"-given conduction status of sw and whether current is" 

" continuous or not, establish derivative of iL" 

PROCEDURAL(piL=sw,contmode, Vcf,Vo, iL,isw) 

!"-if sw gated, assume can conduct pos current" 

IF (sw) THEN 

contmode = .true. 

piL * (-Rlser*iL+Vcf-Vswdrop-Vo)/Lbuck 
ELSE 

"-if iL GT zero, diode conducts" 

IF (contmode) THEN 

piL = (-Riser*iL~Vdidrop-Vo)/Lbuck 
ELSE 

piL = 0.0 !"discontinuous mode" 

ENDIF 

ENDIF 
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END 


! "of procedural" 


!"link filter state variables integrated" 

Vcf = INTEG ( (iLf-isw) /Cf, Vcfic) 

iLf - INTEG((Vdc*voltscal-Vcf-(Rsrc*iLf))/Lf, iLfic) 

!"establish the derivative of Buck output voltage" 
pVo = (iL - io)/Cbuck 

{"integrate the state variables" 

iLub = INTEG(piL, iLic) 

iL = BOUND(-1.Oe-6, 1.0e6, iLub) 

io = Vo/Rload 

Vo = INTEG(pVo, Voic) 

"-form the reference ramp waveform (0 to 1)" 

tt = mod(t, Tramp) 
reframp = tt/Tramp 

"-schedule discontinuous mode when iL tries go neg" 

SCHEDULE DCM .XN. iL 


!"Control Algorithm for Buck" 
CONSTANT hv = 0.00914 
{CONSTANT hi - -0.003299 
CONSTANT hi = -0.003208 
CONSTANT hn = 6.093 


CONSTANT nic =0.0 
CONSTANT intbnd = 1000.0 
CONSTANT minduty =0.0 
CONSTANT maxduty = 0.99 
CONSTANT Kscv = 1.000 
CONSTANT Ksci = 1.000 
va = Vref-Vo 
ioave = Ksci*io 
iLave = Ksci*iL 
ia = iLave-ioave 


! "integral of volt err" 


{"scale factor for voltage sensor" 
J"scale factor for current sensor" 
{"proportional voltage term" 


{"proportional current term" 


pn = va 

nub = INTEG(pn, nic) {"integral voltage term" 

*n = BOUND(-intbnd, intbnd, nub) 
dss = Vref/Vcf 

dub =dss + Kscv*hv*va + Ksci*hi*ia + Kscv*hn*n 
duty = BOUND(minduty, maxduty, dub) 


END {"of derivative" 


vout = Rload*io 

” establish transition to discontinuous mode Buck" 
DISCRETE DCM 

contmode = .false. 

END {"of dcm" 


DISCRETE rampoff 

rampon = .FALSE. 

END 
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END !"of dynamic" 

END !"of program" 

SUBROUTINE simpsonq{thetarsimp,Kone,psiqpripri,psidpripri, iLf 
+ ,iqsrout) 

************************************************************************ 

* 

* subroutine to calculate components for simpson method of determining 

* the q commutating currents 

* 

************************************************************************ 
COMMON / forsub / rad30,rad60,radl20,Xqpripri,Xdpripri,PI, 

+ currscal 

* 

* declare variables 

* 

real ias, parti, part2, part3, Kone, iqsrout, iLf, currscal 

* 

* break the equation for ias down into manageable parts 

* 

partl-Kone-sqrt(3.)*(psiqpripri*cos(thetarsimp+rad30)+psidpripri 
+ *sin(thetarsimp+rad30)) 

* 

part2=(Xqpripri-Xdpripri)*cos(2.*thetarsimp+rad60)+Xqpripri 
+ +Xdpripri 

* 

part3=((Xqpripri-Xdpripri)*cos(2.*thetarsimp+radl20)+0.5 
+ *(Xqpripri+Xdpripri))*(iLf/currscal) 

* 

* combine the 3 parts to determine ias 

* 

ias=(partl/part2)-(part3/part2) 

* 

* now solve for iqsr at the inputed thetar 

* 

iqsrout=( (2.*sqrt(3.))/3.)*(ias*cos(thetarsimp+rad30)- 
+ (iLf/currscal)*sin(thetarsimp)) 

* 

end 

SUBROUTINE simpsond(thetarsimp,Kone,psiqpripri,psidpripri,iLf 
+ ,idsrout) 

************************************************************************ 

* 

* subroutine to calculate components for simpson method of determining 

* the d commutating currents 

* 

************************************************************************ 
COMMON / forsub / rad30,rad60,radl20,Xqpripri,Xdpripri,PI, 

+ currscal 

* 

* declare variables 

* 

real ias, parti, part2, part3, Kone, idsrout, iLf, currscal 

★ 

* break the equation for ias down into manageable parts 

★ 
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partl=Kone-sqrt (3.) * (psiqpripri*cos (thetarsimp+rad30) +psidpripri 
+ *sin (thetarsimp+rad30) ) 

* 

part2=(Xqpripri-Xdpripri)*cos(2.*thetarsimp+rad60)+Xqpripri 
+ +Xdpripri 

* 

part3=((Xqpripri-Xdpripri)*cos(2.*thetarsimp+radl20)+0.5 
+ *(Xqpripri+Xdpripri))*(iLf/currscal) 

* 

* combine the 3 parts to determine ias 

★ 

ias=(partl/part2)-(part3/part2) 

* 

* now solve for idsr at the inputed thetar 

* 

idsrout=((2•*sqrt(3.))/3.)*(ias*sin(thetarsimp+rad30)+ 

+ (iLf/currscal)*cos(thetarsimp)) 


REAL FUNCTION findu(Kone,psiqpripri,psidpripri,iLf) 


function subprogram to calculate u using the ’half-it' method 


COMMON / forsub / rad30,rad60,radl20,Xqpripri,Xdpripri,PI, 

+ currscal 

* 

* declare variables to be used 

* 

logical condition 

real utemp, fcom, fpricom, Kone, iLf, ul, u2, u, currscal 

* 

ul=0.0 
u2=rad60 
u=0.0 

fcom=Kone-sqrt(3.)*(-psiqpripri*sin(beta+u)+psidpripri 
+ *cos(beta+u))-((Xqpripri-Xdpripri)*cos(2.*beta+2.*u 
+ -radl20)+0.5*(Xqpripri+Xdpripri))*(iLf/currscal) 

10 if (ABS(ul-u2) .gt. 0.001) then 
u = 0.5*(ul+u2) 

fcom=Kone-sqrt(3.)*(-psiqpripri*sin(beta+u)+psidpripri 
+ *cos(beta+u))-((Xqpripri-Xdpripri)*cos(2.*beta+2.*u 

+ -radl20)+0.5*(Xqpripri+Xdpripri))*(iLf/currscal) 

if (fcom .le. 0.00) then 
ul=u 
else 
u2=u 
endif 
goto 10 
endif 
findu=u 
END 

REAL FUNCTION findbeta(psidpripri,psiqpripri,iLf, alfahat) 


* function subprogram to calculate beta 
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COMMON / forsub / rad30,rad60,radl20,Xqpripri,Xdpripri,PI, 

+ currscal 

★ 

* declare variables to be used 

* 

real betatop, betabot, betamin, alfahat, fbmin, iLf, currscal 

★ 

betatop = PI 
betabot = -PI 
betamin = 0.0 

* 

* calculate the function to determine if beta is greater than betamin 

* 

fbmin = SQRT(3.)*(psiqpripri*cos(betamin)+psidpripri*sin(betamin 
+ +2. *(iLf/currscal)*(Xqpripri-Xdpripri) 

+ *sin(2.*betamin-rad60) 

10 if (ABS(betatop-betabot) . gt. 0.001) then 
betamin = (betatop+betabot)/2. 

fbmin = SQRT(3.)*(psiqpripri*cos(betamin)tpsidpripri 
+ *sin(betamin))+2.*(iLf/currscal)*(Xqpripri-Xdpripri) 

+ *sin(2.*betamin-rad60) 

if (fbmin .le. 0.0) then 
betabot = betamin 
else 

betatop = betamin 
endif 
goto 10 
endif 

* 

* now determine beta based on the above determination 

* 

findbeta = alfahat + betamin 
END 


B. .CMD FILE 

s strplt = .t. 
s calplt = .f. 

s devplt = 1 

s ppoplt = .f. 
s xinspl = 6 
s weditg = .f. 
s wesitg = .f. 
s nrwitg = .f. 


! "one variable per x-axis" 

! "1 for X-windows" 

! "5 for postscript” 

! "true rotates plot 90 deg" 

! "x-axis plot units" 

! "false suppresses data write" 

! "each time SCHEDULE occurs" 

! "true enables accumulation of data" 



! "after a CONTIN" 


s alcplt = .f. ! "false makes all plots dark" 

s tstop = 1.0 
s maxt=5.0e-7 
s cint=1.0e-3 


prepare 

t,Vdc,Te,psikqpri,psikdpri,psifdpri,iqsr,idsr,u,thetar,beta,alpha,phi& 
Vcf,iLf,duty,Vo,io,exfd 

! "runs a start up to 100%, downpower to 10%, and uppower to 100%" 

PROCED runsystem 

s tstop = 0.5 
start 

s cint = le-4 
s devplt = 5 
s pit = 1 

plot /xlo=0.0 /xhi=0.5 Vdc,Vcf,Vo 
s pit = 2 

plot /xlo=0.0 /xhi=0.5 iLf,io 
s pit = 1 

plot /xlo=0.0 /xhi=0.5 Vdc,Vcf,Vo 
s pit = 2 

plot /xlo=0.0 /xhi=0.5 iLf,io 
s devplt = 1 
s tstop = 0.52 
contin . 

s nrwitg = .t. 
s Rmax = 45.125 
s tstop = 0.57 
contin 

s Rmax = 4.5125 
s tstop =0.62 
contin 

s devplt = 5 
s pit = 4 

plot /xlo=0.5 /xhi=0.62 Vdc,Vcf,Vo 
s pit = 10 

plot /xlo=0.5 /xhi=0.62 iLf,io 
s pit = 4 

plot /xlo=0.5 /xhi=0.62 Vdc,Vcf,Vo 
s pit = 10 

plot /xlo=0.5 /xhi=0.62 iLf,io 
s devplt = 1 
s nrwitg = .f. 
s tstop =0.65 
contin 

s tstop =0.67 
s Cf = 400e-6 
contin 

s nrwitg = .t. 
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s Rmax = 45.125 
s tstop = 0.72 
contin 

s Rmax = 4.5125 


s tstop = 0.77 
contin 
s pit = 13 
s devplt = 5 
plot /xlo=0.65 
s pit = 14 
plot /xlo=0.65 
s pit = 13 
plot /xlo=0.65 
s pit = 14 
plot /xlo=0.65 

END 


/xhi=0.77 Vdc,Vcf,Vo 
/xhi=0.77 iLf,io 
/xhi-0.77 Vdc,Vcf,Vo 
/xhi=0.77 iLf,io 


!"runs a start up to 10%, uppower to 100% and downpower to 10%" 

PROCED runsystem2 

s tstop = 0.5 
s Rmax = 45.125 
start 

s cint = le-4 
s devplt = 5 
s pit = 14 

plot /xlo=0.0 /xhi=0.5 Vdc,Vcf,Vo 
s pit = 15 

plot /xlo=0.0 /xhi=0.5 iLf,io 
s pit = 14 

plot /xlo=0.0 /xhi=0.5 Vdc,Vcf,Vo 
s pit = 15 

plot /xlo=0.0 /xhi=0.5 iLf,io 
s devplt = 1 
s tstop = 0.52 
contin 

s nrwitg = .t. 
s Rmax = 4.5125 
s tstop = 0.57 
contin 

s Rmax = 45.125 
s tstop = 0.62 
contin 

s devplt = 5 
s pit = 16 

plot /xlo=0.5 /xhi=0.62 Vdc,Vcf,Vo 
s pit = 17 

plot /xlo=0.5 /xhi=0.62 iLf,io 
s pit =16 

plot /xlo=0.5 /xhi=0.62 Vdc,Vcf,Vo 
s pit = 17 

plot /xlo=0.5 /xhi=0.62 iLf,io 
s devplt = 1 
s nrwitg = .f. 
s tstop = 0.65 
contin 
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s tstop = 0.67 
s Cf = 4 00e-6 
contin 

s nrwitg = .t. 
s Rmax = 4.5125 
s tstop = 0.72 
contin 

s Rmax =45.125 
s tstop = 0.77 
contin 
s pit = 18 
s devplt = 5 

plot /xlo=0.65 / xhi=0.77 Vdc^cf^Vo 
s pit = 19 

plot /xlo=0.65 /xhi=0.77 iLf,io 
s pit = 18 

plot /xlo=0.65 /xhi=0.77 Vdc f Vcf,Vo 
s pit = 19 

plot /xlo=0.65 /xhi==0.77 iLf,io 


END 
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